/* Homework 2, problem 1.  Compile with
 *
 *	make hw2_1
 *
 * Run with
 *
 *	hw2_1  N
 *
 * where N is the number of nodes in the network.
 */

/*------------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nrutil.h"
#include "nr.h"


#define N	5			/* default */

/* The system consists of N fully connected nodes.  The redundancy in
 * Kirchhoff's laws means that we need only solve the (N-1)x(N-1) system
 * obtained from consideration of nodes 2 through N.  These will be
 * labeled by the indices i and j running from 1 to N-1 below.  (We
 * could also have used dmatrix(2, N, 2, N) below and run from 2 to N.)
 */

double resistance(int i, int j)
{
    if (i > j) {

	/* Replace i and j by min(i,j) and max(i,j), respectively. */

	int tmp = i;
	i = j;
	j = tmp;
    }
    return i + 2*j;
}

int main(int argc, char *argv[])
{
    int i, j, k;
    double **a, **b;
    int nres = N, n1;

    if (argc > 1) nres = atoi(argv[1]);
    n1 = nres - 1;

    /* NR declarations of matrices: */

    a  = dmatrix(1, n1, 1, n1);
    b  = dmatrix(1, n1, 1, 1);		/* b is really a column vector,
					   but gaussj expects a matrix
					   since we can solve multiple
					   equations with a single call */

    /*--------------------------------------------------------------------*/

    /* Set up A and b. */

    for (i = 1; i <= n1; i++)	     /* i = row,    corresponds to node i+1 */
	for (j = 1; j <= n1; j++) {  /* j = column, corresponds to node j+1 */
	    if (i == j) {
		a[i][i] = 0;
		for (k = 1; k <= nres; k++)
		    if (k != i+1) a[i][i] -= 1/resistance(i+1, k);
	    } else
		a[i][j] = 1/resistance(i+1, j+1);
	}

    b[1][1] = 1;		    /* unit external current, so the desired
				       resistance is V2 - V1 = V2 here. */

    for (i = 2; i <= n1; i++) b[i][1] = 0;

    /*--------------------------------------------------------------------*/

    /* Solve by Gauss-Jordan.  On return, A is replaced by
       its inverse and b is replaced by the solution matrix. */

    gaussj(a, n1, b, 1);

    /*--------------------------------------------------------------------*/

    /* Print out the solution. */

    printf("N = %d, R = -V2 = %f\n", nres, -b[1][1]);

    /*--------------------------------------------------------------------*/

}