#define vstring	"IND version 1.2"
/*
 * Compute interferometer coupling taking into account film thickness,
 * which may be less than the penetration depth.
 * Written by Stephen R. Whiteley pre-1990.
 * 8/28/90 Altered expression for flux under control line, see comment.
 * Updated to avoid compiler warnings 11/7/2014.
 *
 * To build:  cc -o ind ind.c -lm
 *
 * This is unrestricted freeware.
 */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

/* comment this out if compiling for other than an IBM PC or compatible
 * with VGA graphics
#define VGA
#include <dos.h>
 */


#define MU 1.256637

struct input {
	double w1;
	double w2;
	double h1;
	double h2;
	double t0;
	double t1;
	double t2;
	double l0;
	double l1;
	double l2;
} in;

struct output {
	double b1;
	double B1;
	double B2;
	double p1;
	double p2;
	double l1;
	double l2;
	double m;
	double k12;
	double k21;
} out ;

#ifdef __STDC__
extern void getinput(void);
extern void solve(void);
extern void printout(void);
extern void printgeom(void);
extern void vgamode(int);
extern void lubksb(double**,int,int*,double*);
extern void ludcmp(double**,int,int*);
extern void nrerror(char*);
#else
extern void getinput();
extern void solve();
extern void printout();
extern void printgeom();
extern void vgamode();
extern void lubksb();
extern void ludcmp();
extern void nrerror();
#endif

int
main()
{
	char s[8];

#if (__TURBOC__)
	_control87(0xffff,0x003f);
#endif
	in.w1 = 10;
	in.w2 = 5;
	in.h1 = .11;
	in.h2 = .5;
	in.t0 = .13;
	in.t1 = .3;
	in.t2 = .6;
	in.l0 = .085;
	in.l1 = .085;
	in.l2 = .085;

	vgamode(0x12);
	printf("\n%s\n\n",vstring);
	printf("Interferometer coupling model\n\n"); 
	printgeom();
	printf("\n\nEnter geometrical data in microns:\n\n");
	while (1) {
		getinput();
		vgamode(0x11);
		solve();
		printout();
		printf("\nQuit ? [y] ");
		(void) fgets(s,7,stdin);
		if (*s == 'y' || *s == 'Y' || *s == '\n') {
			vgamode(0x14);
			exit(0);
		}
	}
}


void
getinput()

/* prompt for new geometrical parameters */
{
	char s[32];
	double d, atof();

	while (1) {
		*s = '\0';
		printf("enter W1  [%8.4f]: ",in.w1);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.w1;
		if (d > 0 && d < 50) {
			in.w1 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter W2  [%8.4f]: ",in.w2);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.w2;
		if (d > 0 && d < in.w1) {
			in.w2 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter H1  [%8.4f]: ",in.h1);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.h1;
		if (d > 0 && d < 2) {
			in.h1 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter H2  [%8.4f]: ",in.h2);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.h2;
		if (d > 0 && d < 2) {
			in.h2 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter T0  [%8.4f]: ",in.t0);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.t0;
		if (d > 0 && d < 2) {
			in.t0 = d;
			break;
		}
		printf("ERROR re");
	}
	while (1) {
		*s = '\0';
		printf("enter T1  [%8.4f]: ",in.t1);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.t1;
		if (d > 0 && d < 2) {
			in.t1 = d;
			break;
		}
		printf("ERROR re");
	}
	while (1) {
		*s = '\0';
		printf("enter T2  [%8.4f]: ",in.t2);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.t2;
		if (d > 0 && d < 2) {
			in.t2 = d;
			break;
		}
		printf("ERROR re");
	}
	while (1) {
		*s = '\0';
		printf("enter L0  [%8.4f]: ",in.l0);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.l0;
		if (d > 0 && d < 2) {
			in.l0 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter L1  [%8.4f]: ",in.l1);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.l1;
		if (d > 0 && d < 2) {
			in.l1 = d;
			break;
		}
		printf("ERROR re");
	}

	while (1) {
		*s = '\0';
		printf("enter L2  [%8.4f]: ",in.l2);
		(void) fgets(s,32,stdin);
        if (sscanf(s, "%lf", &d) != 1)
            d = in.l2;
		if (d > 0 && d < 2) {
			in.l2 = d;
			break;
		}
		printf("ERROR re");
	}

}


void
solve()

/* solve the equation set */
{
	double **a, *b;
	int indx[3], i;
	double l0, l2, x0, x1, x2;

	l0 = in.h1 + in.l0*cosh(in.t0/in.l0)/sinh(in.t0/in.l0) +
		in.l1*cosh(in.t1/in.l1)/sinh(in.t1/in.l1);

	l2 = in.h2 + in.l1*cosh(in.t1/in.l1)/sinh(in.t1/in.l1) +
		in.l2*cosh(in.t2/in.l2)/sinh(in.t2/in.l2);

	x0 = in.l0/sinh(in.t0/in.l0);
	x1 = in.l1/sinh(in.t1/in.l1);
	x2 = in.l2/sinh(in.t2/in.l2);

	a = (double **) malloc(3*sizeof(double *));
	for (i = 0; i < 3; i++) {
		*(a+i) = (double *) malloc(3*sizeof(double));
		memset(*(a+i),0,3*sizeof(double));
	}

	b = (double *) &out;
	memset(b,0,3*sizeof(double));
	memset(indx,0,3*sizeof(int));

	b[0] = 1.0;

	/* b = b1, B1, B2 */
	/* x = I, 0, 0 */
	/* a*b = x */

	/* current in control line */
	/* I = B2*W2/MU	*/
	a[0][2] = in.w2/MU;

	/* current in second line */
	/* b1*(W1-W2) + B1*W2 - B2*W2 = 0 */
	a[1][0] = in.w1-in.w2;
	a[1][1] = in.w2;
	a[1][2] = -in.w2;

	/* flux conservation under second line */
	/* b1*(l0-x0) - B1*(l0-x0) + B2*x1 = 0 */
	a[2][0] = l0 - x0;
	a[2][1] = x0 - l0;
	a[2][2] = x1;

	ludcmp(a,3,indx);
	lubksb(a,3,indx,b);

	/* Comment :
	 * It is not really clear which path is the appropriate	choice for
	 * evaluating the flux operator.  The choice for self inductance
	 * requires the outer surfaces be neglected to reproduce the
	 * textbook inductance expression.  The mutual inductance (p1)
	 * expression provided by taking two surfaces of the ground plane
	 * and the lower surface of the middle electrode gives an expression
	 * which looks good physically ahd has correct behavior in limits,
	 * unlike expressions derived from other paths.  The p2 path takes
	 * the two control line surfaces plus the upper surface	of the
	 * middle electrode (version 1.2) or both surfaces of the middle
	 * electrode and the lower surface of the control line (version 1.1).
	 * It is not clear which is "correct."
	 */

	out.p1 = out.b1*(l0 - x0);

	/* version 1.1 
	out.p2 = out.B2*(l2 - x1) - out.B1*x1;
	*/
	out.p2 = out.B2*(l2 - x2) - out.B1*x1;

	/* upper ground plane and lower middle */
	out.l1 = MU*l0/in.w1;

	/* lower control, both middle, upper ground plane */
	out.l2 = out.B1*(l0-x1) + out.B2*(l2 - x1);

	out.m  = out.p1;
	out.k21 = out.m/out.l1;
	out.k12 = out.m/out.l2;

	for (i = 0; i < 3; i++)
		free(*(a+i));
	free(a);
}


void
printout()

/* print stuff */
{
	printf("Control line (top) driven by 1 mA current\n\n");
	printf("B field amplitude   Flux in region\n");
	printf("b1 = %10.7f\n",out.b1);
	printf("B1 = %10.7f     P1 = %10.7f\n",out.B1,out.p1);
	printf("B2 = %10.7f     P2 = %10.7f\n",out.B2,out.p2);

	printf("\n");
	printf("Interferometer self inductance      = % 10.7f pH per micron\n",
		out.l1);
	printf("(Control line open)\n\n");

	printf("Effective interferometer inductance = % 10.7f pH per micron\n",
		out.l1-out.m*out.m/out.l2);
	printf("(Control line shorted)\n\n");

	printf("Control line self inductance        = % 10.7f pH per micron\n",
		out.l2);
	printf("(Interferometer line open)\n\n");

	printf("Effective control line inductance   = % 10.7f pH per micron\n",
		out.l2-out.m*out.m/out.l1);
	printf("(Interferometer line shorted)\n\n");

	printf("Mutual inductance                   = % 10.7f pH per micron\n\n",
		out.m);

	printf("K parameter                         = % 10.7f\n\n",
		sqrt(out.m*out.m/(out.l1*out.l2)));
	printf("Current transfer ratio, control line to interferometer = % 10.7f\n",
		out.k21);
	printf("Current transfer ratio, interferometer to control line = % 10.7f\n",
		out.k12);
}


void
printgeom()

/* print a diagram of the geometry */
{

printf("\
		        INTERFEROMETER GEOMETRY\n\n");
printf("\
Penetration:\n");
printf("\
		      |<-          W2            ->|\n");
printf("\
Top        L2\n");
printf("\
Second     L1      __  ____________________________\n");
printf("\
Grnd Plane L0	   ^  |				   |\n");
printf("\
		      |	  			   |\n");
printf("\
		   T2 |	  	    .	      	   |\n");
printf("\
	           __ |____________________________|\n");
printf("\
	           ^\n\
	           H2	 	   B2 ->\n");
printf("\
       __  ______________________________________________________\n");
printf("\
       ^  |		  	    x 				 |\n");
printf("\
	  |            						 |\n");
printf("\
       T1 |   .	                    . 	     		    .	 |\n");
printf("\
       __ |______________________________________________________|\n");
printf("\
       ^\n\
\n\
       H1    b1	->    |            B1 ->            |      b1 ->\n");
printf("\
____________________________________________________________________________\n");
printf("\
       ^      x	          	    x	      		    x\n");
printf("\
       T0\n\
\n\
____________________________________________________________________________\n");
printf("\
\n\
	  |<-	                    W1 			       ->|\n");
}


void
vgamode(mode)

/* set up screen for 25, 28, or 50 lines (VGA only) */
int mode;
{
	/* mode == 0x14 25 lines, 0x12 50 lines, 0x11 28 lines	*/

#ifndef VGA
	return;
#else

	_AL = 3;
	_AH = 0;
	geninterrupt(0x10);	/* set mode		*/

	_AH = 0x11;
	_AL = mode;
	_BL = 0;
	geninterrupt(0x10);	/* load alt char set	*/

	_AH = 0x12;
	_BL = 0x20;
	geninterrupt(0x10);	/* load new prnt_scrn	*/
#endif
}

#define TINY 1e-20

void
lubksb(a,n,indx,b)

/* solve ax=b
 * from "Numerical Recipes in C"
 * a is matrix in LU form, n is size of a, indx points to row swap codes
 * b is input/output vector
 */
double **a;
int n, *indx;
double *b;
{
	int i, ii = -1, ip, j;
	double sum;

	for (i = 0; i < n; i++) {
		ip = indx[i];
		sum = b[ip];
		b[ip] = b[i];
		if (ii != -1)
			for (j = ii; j < i; j++)
				sum -= a[i][j]*b[j];
		else if (fabs(sum) > TINY/2.0)
			ii = i;
		b[i] = sum;
	}
	for (i = n-1; i >= 0; i--) {
		sum = b[i];
		for (j = i+1; j < n; j++)
			sum -= a[i][j]*b[j];
		b[i] = sum/a[i][i];
	}
}


void
ludcmp(a,n,indx)

/* set up LU matrix in place.
 * from "Numerical Recipes in C"
 * a is matrix, n is size of a, indx points to row swap codes
 */
double **a;
int n, *indx;
{
	int i, imax, j, k;
	double big, dum, sum, temp;
	double *vv;

	vv = (double *) malloc(n*sizeof(double));
	for (i = 0; i < n; i++) {
		big = 0.0;
		for (j = 0; j < n; j++)
			if ((temp = fabs(a[i][j])) > big)
				big = temp;
		if (big == 0.0)
			nrerror("Singular matrix in LUDCMP");
		vv[i] = 1.0/big;
	}
	for (j = 0; j < n; j++) {
		for (i = 0; i < j; i++) {
			sum = a[i][j];
			for (k = 0; k < i; k++)
				sum -= a[i][k]*a[k][j];
			a[i][j] = sum;
		}
		big = 0.0;
		for (i = j; i < n; i++) {
			sum = a[i][j];
			for (k = 0; k < j; k++)
				sum -= a[i][k]*a[k][j];
			a[i][j] = sum;
			if ((dum = vv[i]*fabs(sum)) >= big) {
				big = dum;
				imax = i;
			}
		}
		if (j != imax) {
			for (k = 0; k < n; k++) {
				dum = a[imax][k];
				a[imax][k] = a[j][k];
				a[j][k] = dum;
			}
			vv[imax] = vv[j];
		}
		indx[j] = imax;
		if (fabs(a[j][j]) < TINY) {
			if (a[j][j] < 0.0) a[j][j] = -TINY;
			else a[j][j] = TINY;
		}
		if (j != n) {
			dum = 1.0/a[j][j];
			for (i = j+1; i < n; i++)
				a[i][j] *= dum;
		}
	}
	free(vv);
}


void
nrerror(errstr)

/* print error message and die	*/
char *errstr;
{
	vgamode(0x14);
	(void) fprintf(stderr,"%s\n",errstr);
	exit(1);
}