/*
 * SLINE - superconducting microstripline model
 * Stephen R. Whiteley (stevew@srware.com)
 * Implements Chang's algorithm from IBM, with a typo-fix in one of
 * the equations.
 * First public release of very ancient program
 * Version 1.0 6/30/96
 * Updated to avoid compiler warnings 11/7/2014
 * 
 * To build: cc -o sline sline.c -lm
 *
 * This is unrestricted freeware.
 */

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

struct params {
    double lwidth;
    double dthick;
    double lthick;
    double gpthick;
    double dielcon;
    double ldepth;
    double gpdepth; };

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

#define EP 8.85416e-6
#define MU 1.256637
#ifndef PI
#define PI 3.1415926
#endif
#define atanh(x) log(((x)+1)/(1-(x)))/2
#define MAX(x,y) (((x)>(y)) ? (x) : (y))

char *STR1 = "\nVary which input? Enter one of\n\
w   (line width)\n\
d   (dielectric thickness)\n\
t   (line thickness)\n\
h   (ground plane thickness)\n";

char *STR2 = "k   (both line and ground plane thickness)\n\
e   (dielectric constant)\n\
l   (line penetration depth)\n\
g   (ground plane penetration depth)\n\
p   (penetration depth, both line and ground plane)\n";

extern void sline();
extern void help();
extern int  prompt();
extern void sens();

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  = 5;
    tline.dthick  = .5;
    tline.lthick  = .6;
    tline.gpthick = .13;
    tline.dielcon = 4.2;
    tline.ldepth  = .085;
    tline.gpdepth = .085;
    printf("\nSRWare Superconducting Microstripline Model\n\n\
Dimensions are in microns.\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 dielectric thickness [%g]: ",&tline.dthick))
            continue;

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

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

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

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

        if (prompt(1,s,"Enter ground plane penetration depth [%g]: ",&tline.gpdepth))
            continue;

        sline(&tline,&out);
        sens(&tline,&out);

        printf("\nZ = %15g ohms\n",out.Z);
        printf("L = %15.5e pH/um\n",out.L);
        printf("C = %15.5e pF/um\n",out.C);
        printf("T = %15.5e psec/um\n",out.T);

        prompt(2,s,"\nEnter q to quit, o to optimize : ",&delta);
        if (strchr(s,'o')) {
            prompt(2,s,"\nWhich output [z,l,c] ? ",&delta);
            switch (*s) {
            case 'l':
                ov = &out.L;
                c = "L value (pH/um) ? ";
                break;
            case 'c':
                ov = &out.C;
                c = "C value (pF/um) ? ";
                break;
            case 'z':
            default:
                ov = &out.Z;
                c = "Z value (ohms) ? ";
            }
            prompt(2,s,c,&value);
            printf("%s", STR1);
            printf("%s", STR2);
            prompt(2,s,"\n? ",&delta);
            switch (*s) {
            case 'p':
                if (tline.gpdepth != tline.ldepth) {
                    printf("Warning!\nSetting ground plane \
penetration depth to equal the line penetration depth.\n");
                    tline.gpdepth = tline.ldepth;
                }
                iv = &tline.ldepth;
                c = "line and ground plane penetration depth";
                break;
            case 'g':
                iv = &tline.gpdepth;
                c = "ground plane penetration depth";
                break;
            case 'l':
                iv = &tline.ldepth;
                c = "line penetration depth";
                break;
            case 'e':
                iv = &tline.dielcon;
                c = "dielectric constant";
                break;
            case 'k':
                if (tline.gpthick != tline.lthick) {
                    printf("Warning!\nsetting ground plane \
thickness to equal the line penetration depth.\n");
                    tline.gpthick = tline.lthick;
                }
                iv = &tline.lthick;
                c = "line and ground plane thickness";
                break;
            case 'h':
                iv = &tline.gpthick;
                c = "ground plane thickness";
                break;
            case 't':
                iv = &tline.lthick;
                c = "line thickness";
                break;
            case 'd':
                iv = &tline.dthick;
                c = "dielectric thickness";
                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;
                if (*s == 'p')
                    tline.gpdepth = tline.ldepth;
                if (*s == 'k')
                    tline.gpthick = tline.lthick;
                sline(&tline,&out);
                if (fabs((*ov - value)/value) < 1e-6) break;
                if (i > 50) {
                    printf("** NO CONVERGENCE **\n");
                    break;
                }
            }
            printf("\n%s = %g\n",c,*iv);
            printf("\nZ = %15g ohms\n",out.Z);
            printf("L = %15.5e pH/um\n",out.L);
            printf("C = %15.5e pF/um\n",out.C);
            printf("T = %15.5e psec/um\n",out.T);
        }
    }
    return (0);
}

void
sline(tl, out)
struct params *tl;
struct output *out;
{
    double p, pp, eta, lnra, rb, rb0;
    double cap, zz, gind, kappa, aa;

    p    = 1 + tl->lthick/tl->dthick;
    p    = 2*p*p - 1;
    p   += sqrt(p*p - 1);
    pp   = sqrt(p);
    aa   = (p+1)/(2*pp);

    eta = aa * (1 + log(4/(p-1))) - 2*atanh(1/pp);
    eta = pp * (PI*tl->lwidth/(2*tl->dthick) + eta);

    /*
    rb0  = eta + aa * log(MAX(p,eta));
    */
    rb0  = eta + .5*(p+1) * log(MAX(p,eta));

    if (tl->lwidth/tl->dthick >= 5)
        rb = rb0;
    else {
        rb  = -sqrt((rb0-1)*(rb0-p));
        rb += (p+1)* atanh(sqrt((rb0-p)/(rb0-1)));
        rb -= 2*pp * atanh(sqrt((rb0-p)/p/(rb0-1)));
        rb += rb0 + pp*(PI*tl->lwidth)/2/tl->dthick;
        if (tl->lwidth/tl->dthick < 1) help(3);
    }

    lnra    = -2*aa * atanh(1/pp);
    lnra   -= 1 + (PI*tl->lwidth)/2/tl->dthick + log((p-1)/4/p);

    kappa   = 2*tl->dthick * (log(2*rb) - lnra)/tl->lwidth/PI;

    cap   = kappa*tl->lwidth*EP*tl->dielcon / tl->dthick;

    gind  = tl->ldepth * (1/tanh(tl->lthick/tl->ldepth) +
        2*pp/rb /sinh(tl->lthick/tl->ldepth));
    gind += tl->gpdepth / tanh(tl->gpthick/tl->gpdepth);
    gind  = MU/tl->lwidth/kappa * (tl->dthick + gind);

    zz    = sqrt(gind / cap);

    out->K = kappa;
    out->L = gind;
    out->C = cap;
    out->Z = zz;
    out->T = sqrt(gind*cap);
}

void
help(i)
int i;
{
    switch (i) {
    case 1:
        printf("\
Enter the line dimensions in microns, or the relative dielectric constant.\n\
Entering r will repeat prompting, and no input will assign the default value.\n");
        printf("\
The sensitivity analysis numbers to be presented represent the ratio of\n\
output variable percentage change to input variable percentage change");
        return;
    case 2:
        printf("\
The optimization algorithm finds the value of the specified input parameter\n\
to yield the value of the output parameter (Z, C, L, or T) given in response\n\
to the prompt.\n");
        return;
    case 3:
        printf("Warning - this geometry may produce inaccuracy\n");
        return;
    }
}

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
sens(tl,out)
struct params *tl;
struct output *out;
{
    double v, ls, cs, zs, ts;
    struct output oplus, ominus;

    printf("\nSensitivity Analysis\n\
%10s %10s   %10s %10s %10s %10s\n"," ","Value","L sens","C sens","Z sens","T sens");

    v = tl->lwidth;
    tl->lwidth = v + .001*v;
    sline(tl,&oplus);
    tl->lwidth = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("line width %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->lwidth = v;

    v = tl->dthick;
    tl->dthick = v + .001*v;
    sline(tl,&oplus);
    tl->dthick = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("diel thick %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->dthick = v;

    v = tl->lthick;
    tl->lthick = v + .001*v;
    sline(tl,&oplus);
    tl->lthick = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("line thick %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->lthick = v;

    v = tl->gpthick;
    tl->gpthick = v + .001*v;
    sline(tl,&oplus);
    tl->gpthick = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("grnd thick %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->gpthick = v;

    v = tl->dielcon;
    tl->dielcon = v + .001*v;
    sline(tl,&oplus);
    tl->dielcon = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("diel const %10f %s %10f %10f %10f %10f\n",v,"  ",ls,cs,zs,ts);
    tl->dielcon = v;

    v = tl->ldepth;
    tl->ldepth = v + .001*v;
    sline(tl,&oplus);
    tl->ldepth = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("line penet %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->ldepth = v;

    v = tl->gpdepth;
    tl->gpdepth = v + .001*v;
    sline(tl,&oplus);
    tl->gpdepth = v - .001*v;
    sline(tl,&ominus);
    ls = 500*(oplus.L - ominus.L)/out->L;
    cs = 500*(oplus.C - ominus.C)/out->C;
    zs = 500*(oplus.Z - ominus.Z)/out->Z;
    ts = 500*(oplus.T - ominus.T)/out->T;
    printf("grnd penet %10f %s %10f %10f %10f %10f\n",v,"um",ls,cs,zs,ts);
    tl->gpdepth = v;
}