/*
 * Compute the transmission line parameters of a coplanar waveguide at
 * an air/substrate interface.  Written by Stephen R.  Whiteley
 * decades ago.  This implements a handbook formula of unknown origin. 
 * Code updated to avoid compiler warnings 11/7/2014.
 *
 *  To build: cc -o cpw cpw.c -lm
 *
 * This is unrestricted freeware.
 */

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


struct params {
	double lwidth;
	double space;
	double dthick;
	double lthick;
	double dielcon; };

struct output {
	double L;
	double C;
	double Z;
	double T; };

#define EP 8.85416e-6
#define MU 1.256637
#ifndef PI
#define PI 3.1415926
#endif
int noprint;

extern void cpw();
extern double ellip();
extern int  prompt();
extern void help();

int
main()
{
	struct params tline;
	struct output out;
	double atof(), fabs(), lasti, lasto, *iv, *ov, delta, value, tiv;
	int i;
	char s[32], *c;

#if (__TURBOC__)
	_control87(0xffff,0x003f);
#endif
	tline.lwidth  = 300;
	tline.space   = 35;
	tline.dthick  = 254;
	tline.lthick  = 1;
	tline.dielcon = 4.2;
	printf("\nCoplanar Transmission Line Model\n\nValues are \
in microns, picofarads, picohenries, ohms, picoseconds.\n");
	printf("Enter q anytime to quit, h for help.\n");
	while (1) {
		if (prompt(1,s,"\nEnter line width [%g]: ",&tline.lwidth))
			continue;

		if (prompt(1,s,"Enter gap width [%g] : ",&tline.space))
			continue;

		if (prompt(1,s,"Enter dielectric thickness [%g]: ",&tline.dthick))
			continue;

		if (prompt(1,s,"Enter line thickness [%g]: ",&tline.lthick))
			continue;

		if (prompt(1,s,"Enter relative dielectric constant [%g]: ",&tline.dielcon))
			continue;

		cpw(&tline,&out);

		if (!noprint) {
			printf("\nZ = %g\n",out.Z);
			printf("C = %e\n",out.C);
			printf("L = %e\n",out.L);
			printf("T = %e\n",out.T);
		}

		prompt(2,s,"\nEnter q to quit, o to optimize : ",&delta);
		if (strchr(s,'o')) {
			prompt(2,s,"which output [z,l,c] ? ",&delta);
			switch (*s) {
			case 'l':
				ov = &out.L;
				c = "L value ? ";
				break;
			case 'c':
				ov = &out.C;
				c = "C value ? ";
				break;
			case 'z':
			default:
				ov = &out.Z;
				c = "Z value ? ";
			}
			prompt(2,s,c,&value);
			prompt(2,s,"vary which input [w,s,e] ? ",&delta);
			switch (*s) {
			case 's':
				iv = &tline.space;
				c = "gap width";
				break;
			case 'e':
				iv = &tline.dielcon;
				c = "dielectric constant";
				break;
			case 'w':
			default:
				iv = &tline.lwidth;
				c = "line width";
			}
			i = 0;
			lasti = 0;
			lasto = 0;
			while (++i) {
				delta = (*iv - lasti)/(*ov - lasto);
				lasto = *ov;
				lasti = *iv;
				tiv   = (value - *ov)*delta;
				if (fabs(tiv) > 0.2*fabs(*iv))
					tiv = ((tiv > 0) ? 0.2*fabs(*iv) :
						-0.2*fabs(*iv));
				*iv += tiv;
				cpw(&tline,&out);
				if (fabs((*ov - value)/value) < .0001) break;
				if (i > 50) {
					printf("** NO CONVERGENCE **\n");
					break;
				}
			}
			if (!noprint) {
				printf("\n%s = %g\n",c,*iv);
				printf("Z = %g\n",out.Z);
				printf("C = %e\n",out.C);
				printf("L = %e\n",out.L);
				printf("T = %e\n",out.T);
			}
		}
	}
}

void
cpw(tl,out)
struct params *tl;
struct output *out;
{
	double k, ke, f, fe, a1, a2, epre, ellip();

	ke  = (1.25*tl->lthick/PI) * (1+log(4*PI*tl->lwidth/tl->lthick));
	ke  = (tl->lwidth + ke) / (2*tl->space + tl->lwidth - ke);

	k   = tl->lwidth/(tl->lwidth+2*tl->space);

	noprint = 0;
	if (ke >= 1) {
		printf("These parameters are not calculable\n");
		noprint = 1;
		return;
	}
	f = ellip(k*k)/ellip(1-k*k);
	fe = ellip(ke*ke)/ellip(1-ke*ke);
	
	a1 = tanh(0.775 * log(tl->dthick/tl->space) + 1.75);
	a2 = 0.04 - 0.7*k + 0.01*(1 - 0.1*tl->dielcon)*(0.25 + k);
	epre = (tl->dielcon+1)* 0.5 *(a1 + k*tl->space*a2/tl->dthick);
	epre = epre - (0.7*(epre-1)*tl->lthick/tl->space) /
		(f + 0.7*tl->lthick/tl->space);

	out->C = 4*EP*epre*fe;
	out->Z = sqrt(MU/EP)/4/sqrt(epre)/fe;
	out->L = out->Z*out->Z*out->C;
	out->T = sqrt(out->L*out->C);
}

double
ellip(y)
double y;
{
	int i, ms;
	double cf, ckk, err, x, fabs();
	x = cf = ckk = 1.0;
	ms = 2000;
	if (fabs(y) > 1) ms = 20;
	for (i = 2; i <= ms; i += 2) {
		cf  *= (double)(i-1)/i;
		ckk *= y;
		err  = cf*cf*ckk;
		x   += err;
		if (fabs(err) <= 1.0e-7) break;
	}
	x *= PI/2;
	return x;
}

int
prompt(i,s,c,v)
int i;
char *s, *c;
double *v;
{
    double d;
    for (;;) {
		printf(c,*v);
        *s = '\0';
        (void)fgets(s, 32, stdin);
		if (strchr(s,'q')) exit(0);
		if (strchr(s,'r')) return 1;
		if (strchr(s,'h')) { help(i); continue; }
        if (sscanf(s, "%lf", &d) == 1)
            *v = d;
		return 0;
	}
}

void
help(i)
int i;
{
	switch (i) {
	case 1:
		printf("\
Enter the line parameters in microns, or the relative dielectric constant.\n\
Entering r will repeat prompting.\n");
		return;
	case 2:
		printf("\
The optimization algorithm finds the value of the specified input parameter\n\
(line width, gap width, or dielectric constant) to yield the value of the\n\
output parameter (Z, C, L, or T) prompted for.\n");
		return;
	case 3:
		printf("Warning - this geometry may produce inaccuracy\n");
		return;
	}
}