/* ------------------------------------------------------ */
/* Maccready theory by backward induction on Mcready value*/
/* This program implements all calculations for           */
/* "MacCready theory with uncertain lift and limited      */
/*  altitude"                                             */
/* The program is a GAUSS program. Every item is a matrix */
/* and dot (.) operators operate element by element.      */
/* The MacCready calculations start on l. 322. To see how */
/* it's done, start there.                                */

/* John H. Cochrane May 4 1999
   john.cochrane@gsb.uchicago.edu                         */
/* ------------------------------------------------------ */

library pgraph; /* initializes graphics routines     */
output file = newmc.out reset;

/* ------------------------------ */
/* switches for how program runs  */
/* -----------------------------  */

justgrf = 0;    /* allows you to skip everything and play with graph at end */
                /* IF you have alreday run once so the results (lamfin) are in
                   the workspace */
thermgf = 0;    /* plot thermal distribution? */
convt   = 0;    /* make rotated files for eps conversion? */
polarplt = 0;   /* plot polar, etc? */
verbose = 0;    /* lots of chat */

                               /* preprogrammed cases */
case = 2;                      /* 1: 0,4 kts, no porpoising, Illinois;
                                                              dry discus  */
                               /* 2: more realistic Illinois; dry discus  */
                               /* 3: 1-26, nimbus, realistic Illinois     */
                               /* 4: wet discus, uvalde                   */
                               /* 0: don't use a preprogrammed case       */

/* -------------------------- */
/* Define glider performance  */
/* -------------------------- */

smin = 1.15;   /* min sink in kts */
vms  = 42;     /* speed for min sink in kts */
vhi = 80;      /* high speed polar point */
shi = 3;       /* high speed sink in kts */
hcap = 0.96;   /* used to scale up or down wining speed assumption */

    /* preprogrammed glider cases. */
    /* 1: dry discus 2: wet Nimbus 3/asw22 3:1-26  4: wet discus*/

if (case == 1) or (case == 2) or (case == 4);
    glider = 1;
endif;
if case == 3;
    glider = 3;
endif;

if glider == 1;
    smin = 1.15;  /* min sink in kts */
    vms  = 42;    /* speed for min sink in kts */
    vhi = 80;      /* high speed polar point */
    shi = 3;      /* high speed sink in kts */
    hcap = 0.96;
endif;
if glider == 2;
    smin = 0.95;  /* min sink in kts */
    vms  = 50;    /* speed for min sink in kts */
    vhi = 100;      /* high speed polar point */
    shi = 3;      /* high speed sink in kts */
    hcap = 0.82;
endif;
if glider == 3;
    smin = 1.7;    /* min sink in kts */
    vms  = 33;    /* speed for min sink in kts */
    vhi = 64;      /* high speed polar point */
    shi = 4;      /* high speed sink in kts */
    hcap = 1.6;
endif;
if glider == 4;
    ldg = (9/6.5)^0.5; /* wingloading factor */
    smin = 1.15*ldg;  /* min sink in kts */
    vms  = 42*ldg;    /* speed for min sink in kts */
    vhi = 80*ldg;      /* high speed polar point */
    shi = 3*ldg;      /* high speed sink in kts */
    hcap = 1;        /* reset V target for this case */
endif;

/* ------------------------------ */
/* Define grids and thermal model */
/* ------------------------------ */

hgrid = seqb(0,5000,500);   /* height grid in feet. MUST be even steps: mixer */
lgrid = 0;                  /* thermal strength grid in kts */
lprb  = 0;                  /* thermal probability in each mile */
porp  = 0;                  /* porpoise factor -- cruse strength/climb */
Xmax  = 150;                /* course length in n mi */
Vwin  = 50;                 /* winner's handicapped speed in kts */
Vwin = Vwin/hcap;
dx    = 1;                 /* step size in n mi */
alpha = 0.65;              /* maximum points for landout */
minclimb = 500;            /* min altitude for thermaling */
fstp = 1;                  /* how many miles back for first step? */
mix = 100;                 /* 1 sigma feet / mile glide uncertainty */


                           /* Preprogrammed cases */
if case == 1;
    hgrid = seqb(0,5000,50);      /* height grid in feet */
    lgrid = (0)|(4+smin);         /* thermal strength grid in kts */
    lprb  = .9|.1;                /* thermal probability in each mile */
    /* use lgrid = 0; lprb = 1; to check that it reverts to glide to finish */
    porp  = 0;                    /* porpoise factor -- cruse strength/climb */
    mix = 50;
endif;

if (case == 2) or (case == 3);
    hgrid = seqb(0,5000,50);
    lgrid = 0|(1+smin)|(2+smin)|(4+smin)|(6+smin);
    lprb  = 0|0.2|0.1|0.05|0.02;
    lprb[1] = 1-sumc(lprb);
    porp  = 0.5;
    mix = 75;
endif;

if case == 4;
    hgrid = seqb(0,9000,50);
    lgrid =  0|   (2+smin)|  (4+smin)|   (6+smin)|(8+smin);
    lprb  =  0|0.15|0.1|0.05|0.025;
    lprb[1] = 1-sumc(lprb);
    porp  = 0.5;
    Vwin = 70;
    Xmax = 250;
endif;

/* ------------------------------ */
/* graphs to assess thermal model */
/* ------------------------------ */

if thermgf;
    graphjc;
    title("Thermal distribution");
    _plctrl = -1;
    _pbarwid = 0.2;
    _pbartyp = 0~15;
    if case == 2;
        _ptek = "thermf2.tkf";
    endif;
    bar(lgrid,lprb);
    if case == 4;
        _ptek = "thermf4.tkf";
    endif;

    if convt;
        barps(lgrid-smin,lprb);
    endif;

    graphjc;
    _pcolor = 15;
    title(
     "Chance finding a thermal as good or better than indicated in x miles");
    prbetr = rev(cumsumc(rev(lprb)));
    pr1 = 1-prbetr;
    px = seqb(1,20,1);
    py = 1-(pr1')^px;
    if (case == 2) ;
        _plegstr = "1 kt\0002 kts\0004 kts\0006 kts";
        _plegctl = 1;
    ENDIF;
    if (case == 4) ;
        _plegstr = "1 kt\0002 kts\0004 kts\0006 kts\0008kts";
        _plegctl = 1;
    ENDIF;
    _pstype = 1|1|2|3;
    _plctrl = 0|1|1|1;
    xtics(0,20,5,0);
    ytics(0,1,0.2,0);
    xlabel("Nautical miles");
    ylabel("Probability");
    if case == 2;
        _ptek = "cumprb2.tkf";
    endif;
    if case == 4;
        _ptek = "cumprb2.tkf";
    endif;
    xy(px,py);
    if convt;
        graphps(px,py);
    endif;
endif;

/* ----------------------------- */
/* Glider performanc procedures  */
/* All procedures use a quadratic*/
/* polar approximation           */
/* ----------------------------- */

/* find coefficient a in s(v) = smin + a/2 (v-vms)^2 to fit s80 */

a = 2*(shi-smin)/(vhi-vms)^2;

/* define Macready and glide functions */

Proc(1) = s(v);
if minc(minc(v)) < 0;
    "s(v) error: asking for speed less than zero v" v;
endif;
retp( (v .>= vms).*(smin+a/2*(v-vms)^2) +
      (v .<  vms).*smin );
endp;

proc(1) = Mc(v);
retp(1/2*a*(v^2 - vms^2) - smin);
endp;

proc(1) = mcinv(lam);
    retp( (lam.>-smin).*(abs(2/a*(lam+smin)+vms^2)^0.5)
          +(lam.<=smin).*0);  /* speed = 0 or climb in this case */
endp;

/* G2 returns inverse glide angle (height lost per mile) */
/* this can handle negatives (porpoising) */
/* speed is restricted to be at least vms so you can't get infinite up */

proc(1) = G2(m,lam);
local vl,gl,laml,ml,i,j;
if (rows(m) > 1) and (rows(lam) > 1);
    "G2 error: one of m, lam may be a vector but not both";
endif;
i = 1;
    do while i <= rows(m);
    ml = m[i];
    j = 1;
    do while j <= rows(lam);
        laml = lam[j];
        vl = mcinv(laml-ml);
        if vl == 0;
            vl = vms;
        endif;
        if (i == 1) and (j == 1);
            gl = (s(vl)-ml)/vl;
        else;
            gl = gl|((s(vl)-ml)/vl);
        endif;
        j = j+1;
    endo;
    i = i+1;
endo;
retp(gl);
endp;

proc(1) = G(m,lam);
local vl,Gl;
vl = mcinv(lam-m);
Gl = vl./(s(vl)-m);
retp(Gl);
endp;

LDmax = G(0,0);

proc(1) = vofg(m,gl); /* returns speed to attain glide g in lift/sink m */
    local al,bl,cl,disc,vl;
    if gl < 0;
        "vofg: being asked for negative glide angle";
    endif;
    if (m .< smin) and (gl > G(m,0));   /* can't achieve this glide angle */
        retp(-1);      /* code for can't get there from here */
    else;
        al = a/2;
        bl = -(1/gl+a*vms);
        cl = smin + a/2*vms^2 - m;
        disc = bl^2 - 4*al*cl;
        if disc < 0 ;
            "vofg error: negative discriminant g l gmax" gl m G(m,0);
        else;
            vl = (-bl + (disc)^0.5)/(2*al);
        endif;
    endif;
    retp(vl);
endp;

/* --------------------------------- */
/* plot polar and related quantities */
/* --------------------------------- */

if polarplt;

graphjc;
v = seqb(30,130,1);
title("sink rate");
xtics(0,120,20,1);
xlabel("speed in kts");
xy(v,s(v));

title("Macready to fly each speed in still air");
xy(v,Mc(v));

title("glide angle");
xy(v,v./s(v));


m = 2|0|-2;
lam = seqb(0.5,10,0.5);
Ga = G(m',lam);

graphjc;
title("Glide as function of MacCready and lift/sink");
_plegstr = "m=+2\000m=0\000m=-2";
_plegctl = 1;
xlabel("McReady");
ytics(-20,60,10,0);
xy(lam,Ga);

graphjc;
title("Speed as function of MacCready and lift/sink");
_plegstr = "m=+2\000m=0\000m=-2";
_plegctl = 1;
xlabel("MacCready");
xy(lam,mcinv(lam-m'));

endif;


if not justgrf; /* allows you to skip calculations and play with graphs */

/* ==================================================================== */
/*                      START MACCREADY CALCULATIONS                    */
/* ==================================================================== */

/* -------------------------------- */
/* Find values one mile out         */
/* simply fly home, no lift or sink */
/* -------------------------------- */

Twin = xmax/vwin;
                                   /* initialize grids */
W = zeros(rows(hgrid),1);          /* speed at finish */
Wh = zeros(rows(hgrid),1);         /* Wh (h,x) */
wt = -ones(rows(hgrid),1)/Twin;    /* Wt (h,x) if no landout */
wt[1] = 0;                         /* at zero, you've landed out */
W[1] = alpha*(Xmax-fstp*dx)/xmax;
                               /* landout rules, assuming winning speedtonow*/
Wh[1] = -999;                  /* code for not known yet */
hi = 2;
do while hi <= rows(hgrid);
    h = hgrid[hi];
    if verbose; "  height" h ; endif;
    v = vofg(0,(6000*fstp*dx)/h);
                         /* speed to just exhaust altitude in this lift*/
    if verbose; "  speed" v ; endif;
    if v >= vms;           /* can make it home */
        if verbose; "  can make it home mc" mc(v); endif;
        W[hi] = xmax/(Twin*(xmax-fstp*dx)/xmax + fstp*dx/v) /Vwin;
                           /* speed at finish, assuming you've flown
                              at winning speed so far */
        wh[hi] = 1/Twin * 1/ ( mc(v) );
    elseif v == -1;        /* can't make it home */
        wt[hi] = 0;
        w[hi] = -999;
        wh[hi] = -999;
    else;
        "error: speed " v " reported to get home";
    endif;
    hi = hi+1;
endo;

if verbose;
    " h  W Wh Wt before interpolation ";
    hgrid[1:10]~W[1:10]~wh[1:10]~Wt[1:10];
endif;

/* now fill in Wh for landout h values by linearly interpolating W */

i = 2;
do while W[i] == -999;
    i = i+1;
endo;
   /* now i is the smallest index of a valid height */
if verbose; "i index" i; endif;

do while i <= rows(hgrid);
    Testwh = 6000*(W[i]-W[1])/(hgrid[i]-hgrid[1]);
                             /* Guess at Wh in landout region */
    if testwh > Wh[i];       /* Don't want a huge peak in Wh */
        break;               /* (infinite value of an inch) */
    else;                    /* hence use linear W until we reach a */
        i = i+1;             /* reasonable Wh from the fly home rule */
    endif;
endo;

wh[1:i-1] = testwh*ones(i-1,1);
W[1:i-1] = W[1]+(hgrid[1:i-1]-hgrid[1])./(hgrid[i]-hgrid[1])*(W[i]-W[1]);
if verbose;
    " h  W Wh Wt before interpolation ";
    hgrid[1:10]~W[1:10]~wh[1:10]~Wt[1:10];
endif;

if verbose;
    graphjc;
    _psymsiz = 1; _plctrl = 1;
    title("W one mile out");
    xy(hgrid,w);
    title("wh one mile out");
    xy(hgrid,wh);
    title("wt one mile out");
    xy(hgrid,wt);
    title("MacCready one mile out");
    py = -wt./wh;
    px = hgrid*ones(1,cols(py));
    _pltype = 6*ones(cols(py),1);
    px = fill(px,g2(0,py)*fstp*6000);
    py = py~py;
    _plegstr = "calculation\000Macready to finish";
    _plegctl = 1;
    _psymsiz = 1;
    _plctrl = 1;
    xy(px,py);
endif;


/* take average over lift values */

whv = wh;
wtv = wt;
lamv = -wtv./whv;
lamfin = lamv*ones(1,fstp); /* stores final answer */

/* we're done! we have Wh and Wt one mile out. */

/* ----------------  */
/* now iterate back  */
/* ----------------- */

whl = zeros(rows(hgrid),rows(lgrid));         /* wh(h,x,l) */
wtl = zeros(rows(hgrid),rows(lgrid));
Whl[1,.] = alpha*ldmax/xmax*ones(1,rows(lgrid)); /* landout value */

holdg = hgrid;                    /* place to store height at x -dx */
xi = fstp+1;                      /* index of distance */
do while xi <= xmax+1;

    /* mix up h */                /* adds subtracts random amount to h */
    whvnew = mixer(whv,hgrid,mix) ;
    wtvnew =  mixer(wtv,hgrid,mix);
    if verbose;
        "Starting mile " Xmax-xi "hit any key to continue"; wait;;
        graphjc;
        title("old and smoothed wh, wt");
        _psymsiz = 1;
        _plctrl = 1;
        xy(hgrid,whv~whvnew~wtv~wtvnew);
    endif;
    whv = whvnew;
    wtv = wtvnew;

    li = 1;                         /* look over all lift values */
    do while li <= rows(lgrid);
        l = lgrid[li];
        hi = 1;                     /* look over all altitudes */
        do while hi <= rows(hgrid);
            h = hgrid[hi];
            /* for each altitude and lift, we have the mcready, so we can work
               back and find the altitude at which you leave to get here */
            lam = lamv[hi];
            holdg[hi] = h + g2(l*porp,lam)*dx*6000;
            hi = hi+1;
        endo;

        /* at this point, whl and wtl are wh and wt, but at the new h.
        next task is to interpolate this to give wh and wt on the old h grid*/

        hi = 2;   /* one always means land out */
        do while hi <= rows(hgrid);
            h = hgrid[hi];
            if h <= minc(holdg); /* you're out */
                whl[hi,li] = alpha*ldmax/xmax;
                wtl[hi,li] = 0;
            else;
                whl[hi,li] = interp(h,holdg,whv);   /* interpolates */
                wtl[hi,li] = interp(h,holdg,wtv);
            endif;
            hi = hi+1;
        endo;

        if verbose;
            graphjc;
            title("1 working back; cruise whl and wtl and interpolation");
            px = fill(holdg*(1~1),hgrid*ones(1,2));
            py = fill(whv~wtv,whl[.,li]~wtl[.,li]);
            _plegstr = "wh before\000wt before" $+
                       "wh after\000wtafter";
            _plegctl = 1;
            _psymsiz = 1;
            _plctrl = 1;
            _pstype = 4;
            xy(px,py);
            title("2 working back; cruise lam and interpolation");
            px = fill(holdg*(1),hgrid*ones(1,1));
            py = fill((-wtv./whv),(-wtl[.,li]./whl[.,li]));
            _plegstr = "lam before\000lam after";
            _plegctl = 1;
            _psymsiz = 1;
            _plctrl = 1;
            _pstype = 4;
            xy(px,py);

        endif;

        /* at this point, we have wht and wlt assuming that you will
           glide. However, you have the option to climb. If l-smin > mc
           then you climb up to the point that l-smin = mc (mc rises with
           altitude) or the top of the lift, whichever comes first. */

         if maxc(whl[.,li] .== 0) ; "error: whl = 0 whl" whl[.,li]'; endif;
         lam = -wtl[.,li]./whl[.,li]; /* MacCready if you criuse */
         flg = 0;               /* flag for, have you foud point to leave yet */
         if minc(lam) < l - smin; /*if there are any altitudes where you climb*/

                                  /* find altitude where you stop thermaling */
            indx = rows(hgrid);   /* start at top */
            do while hgrid[indx] > minclimb; /* can't thermal below minclimb */
                if lam[indx] < l - smin;     /* if lift too weak, leave here */

                    if flg == 0;             /* if this is where you leave */
                        wtall = wtl[indx,li];
                        if indx == rows(hgrid);   /* if leave at top, wh jumps*/
                            whall = -wtall/(l - smin);
                        else;
                            whall = -wtall/(l - smin); /*let wh jump to soak
                                                         approx error */
                        endif;
                        flg = 1;
                    endif;
                   if flg == 1;             /* if you will leave higher, Mc work
                                               down */
                        wtl[indx,li] = wtall;
                        whl[indx,li] = whall;
                    endif;
                endif;
                indx = indx-1;
            endo;
         endif;

        li = li+1;                          /* go on to next lift */
    endo;
    whv = sumc(lprb.*whl');                 /* expected values over lift */
    wtv = sumc(lprb.*wtl');
    lamv = -wtv./whv;
    xi = xi+1;                              /* move back one mile */
    lamfin = lamfin~lamv;                   /* store final results */
endo;

endif; /* ends justgrf */


/* ----------------------------------------------- */
/* Graph MacCready rules vs. height, distance to go  */
/* ----------------------------------------------- */

if case == 1;
   graphjc;
    _pcolor = 15;
    xlabel("MacCready setting");
    ylabel("Altitude");
    _plegstr = "10 nm out\00020\000100\000150";
    _plegctl = 0;

    _Pmsgstr = "10 nm out\00020 nm out\000100\000150";
    _pmsgctl =  (4  ~ 1700 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
                (4  ~ 3800 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
                (3.2  ~ 5000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
                (2.7  ~ 5000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
 /* _Pmsgctl =  x(ll) ~ y(ll) ~ ht(.15)  ~angle(o) ~  plot(1)  ~ co(10)~ thick*/

    xtics(0,5,1,0);
    ytics(0,5000,1000,1);
    px = lamfin[.,10|20|100|150];
    py = hgrid;
    _pltype = 6|6|6|6|3|3|3|3|3;
    py = py*ones(1,cols(px));
    mcvals = seqb(0,5,.1);
    px = fill(px,mcvals*ones(1,2));
    py = fill(py,g2(0,mcvals)*(10~20)*6000);
    _ptek = "case1.tkf";
    xy(px,py);
    if convt;
        graphps(px,py);
    endif;


endif;


if  case == 2;
   graphjc;
    _pcolor = 15;
    xlabel("MacCready setting");
    ylabel("Altitude");
    _plegstr = "10 nm out\00020\00030\000100\000150";
    _plegctl = 0;

 _Pmsgstr = "10 nm out\00020 nm out\000100\000150";
 _pmsgctl =  (4  ~ 1600 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (4  ~ 3300 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (3.8  ~ 4800 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (3.2  ~ 4800 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
@ _pmsgctl = 0;@
/* _Pmsgctl =  x(ll) ~ y(ll) ~ ht(.15)  ~angle(o) ~  plot(1)  ~ co(10)~ thick*/

    xtics(0,6,1,0);
    ytics(0,5000,1000,1);
    px = lamfin[.,10|20|100|150];
    py = hgrid;
    _pltype = 6|6|6|6|3|3|3|3|3;
    py = py*ones(1,cols(px));
    mcvals = seqb(0,6,.2);
    px = fill(px,mcvals*ones(1,2));
    py = fill(py,g2(0,mcvals)*(10~20)*6000);
    _ptek = "case2.tkf";
    xy(px,py);
    if convt;
        graphps(px,py);
    endif;

endif;

if  case == 3;
   graphjc;
    _pcolor = 15;
    xlabel("MacCready setting");
    ylabel("Altitude");
    _plegstr = "10 nm out\00020\00030\000100\000150";
    _plegctl = 0;

 _Pmsgstr = "10 nm out\00020\000100\000150";
 _pmsgctl =  (4  ~ 3000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (1.5  ~ 4600 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (2.8 ~ 4600 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
@ _pmsgctl = 0;@
/* _Pmsgctl =  x(ll) ~ y(ll) ~ ht(.15)  ~angle(o) ~  plot(1)  ~ co(10)~ thick*/

    xtics(0,6,1,0);
    ytics(0,5000,1000,1);
    px = lamfin[.,10|20|100];
    py = hgrid;
    _pltype = 6|6|6|3|3|3|3|3;
    py = py*ones(1,cols(px));
    mcvals = seqb(0,6,.1);
    px = fill(px,mcvals*ones(1,2));
    py = fill(py,g2(0,mcvals)*(10~20)*6000);
    _ptek = "case3.tkf";
    xy(px,py);
    if convt;
        graphps(px,py);
    endif;

endif;

if  case == 4;
   graphjc;
    _pcolor = 15;
    xlabel("MacCready setting");
    ylabel("Altitude");
    _plegstr = "20 nm out\00050\000250";
    _plegctl = 0;

 _Pmsgstr = "20 nm out\00035\00050\000250";
 _pmsgctl =  (6  ~ 3900 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (6.8  ~ 7000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (6.4  ~ 8000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 )|
             (5.5  ~ 8000 ~ 0.2 ~ 0 ~ 1 ~ 15 ~ 1 );
@ _pmsgctl = 0;@
/* _Pmsgctl =  x(ll) ~ y(ll) ~ ht(.15)  ~angle(o) ~  plot(1)  ~ co(10)~ thick*/

    xtics(0,8,1,0);
    ytics(0,9000,1000,1);
    px = lamfin[.,20|35|50|250];
    py = hgrid;
    _pltype = 6|6|6|6|3|3|3|3|3;
    py = py*ones(1,cols(px));
    mcvals = seqb(0,8,.1);
    px = fill(px,mcvals*ones(1,2));
    py = fill(py,g2(0,mcvals)*(20~35)*6000);
    _ptek = "case4.tkf";
    xy(px,py);
    if convt;
        graphps(px,py);
    endif;

endif;

/* ========================================================== */
/*                       PROCEDURES                           */
/* ========================================================== */

/* interpolation proc */
/* assumes x is low to high */

proc(1) = interp(x,xgrid,ygrid);
    local T,indxv,dif,low,high,indxl,indxh ,xl,xh,yl,yh,y;

    if minc(xgrid[2:rows(xgrid)]-xgrid[1:rows(xgrid)-1]) < 0 ;
        "interp: error. xgrid must be organized from low to high";
    endif;

    T = rows(xgrid);
    if x < xgrid[1];
        y = ygrid[1] - (x-xgrid[1])*(ygrid[2]-ygrid[1])/(xgrid[2]-xgrid[1]);
    elseif x > xgrid[rows(xgrid)];
        y = ygrid[T] + (x-xgrid[T])*(ygrid[T]-ygrid[T-1])/(xgrid[T]-xgrid[T-1]);
    else;

        indxv = seqb(1,rows(hgrid),1);
        dif = x-xgrid;
        low =  dif .>= 0;
        indxl = maxindc(indxv.*low);
        indxh = indxl + 1;
        xl = xgrid[indxl];
        yl = ygrid[indxl];
        xh = xgrid[indxh];
        yh = ygrid[indxh];
        y = (yh-yl)/(xh-xl)*(x-xl)+yl;
    endif;
    retp(y);
    endp;

/* proc to mix up altitudes according to normal distribution */

proc(1) = mixer(y,x,mix);
  local yn,delx,indxs,wts,i,test,wl,il;
  if mix == 0;
        yn = y;
  else;
      yn = zeros(rows(y),1);
      delx = x[2]-x[1]; /* assumes same grid spacing for all height! */
      indxs = seqb(-10,10,1);
      wts = exp(-1/2*((indxs*delx)^2)/(mix^2));
      i = 1;
      do while i <= rows(x);
          test =  (i+indxs .> 0) .and (i+indxs .< rows(x));
          wl = selif(wts,test);
          wl = wl/sumc(wl);
          il = selif(indxs,test);
          yn[i] = sumc(wl.*y[i+il]);
          i = i+1;
      endo;
  endif;
  retp(yn);
  endp;

/*                           GRAPHJC.G
customized version of graphset. Use in place of graphset.
must call library pgraph first.
*/

proc(0) = graphjc;

graphset;
_pdate = "";
_pltype = 6;
_plctrl = 0;
_pnum = 1|2;
_pframe = 0|0;
fonts("simplex simgrma");
/*
  _pcolor = 15; makes it all black for exporting
*/
endp;


/* ----------------------------------------------------   */
/*                      FILL.G                            */
/*  proc fill to put series of different lengths on graph */
/* ----------------------------------------------------   */

proc(1) = fill(x,y);
    if rows(x) < rows(y);
        retp( (x|(ones(rows(y)-rows(x),1).*x[rows(x),.]))~y);
    endif;
    if rows(x) > rows(y);
        retp( x~(y|(ones(rows(x)-rows(y),1).*y[rows(y),.])));
    endif;
    if rows(x) == rows(y);
        retp(x~y);
    endif;
endp;

/* ------------------------------------------------------------------------- */
/*                             SEQB.G                                        */
/*
    An alternative to seqa.  specify instead min, max, and a step size.  If
    max is not min + integer number * step size, the next smaller number is
    used.

    Usage:

    x = seqb(min,max,step);

*/
/* ------------------------------------------------------------------------- */
proc seqb(min,max,step);
    local nsteps, x;
     if max < min ; "seqb: asking for max < min"; endif;
    nsteps = trunc((max-min)/step+1.000000000001);
    x = seqa(min,step,nsteps);
    retp(x);
endp;

/*                  graphps.g

  program to produce graphs rotated on top of page,
  ready for importation into scientific wor
  replaces xy(px,py)
  call graphset of graphjc first, and then add titles, etc

*/
/*
example of use:
*/
/*
library pgraph;
graphjc;

px = seqa(1,1,500)/50;
py = sin(px);
title("Title here");
graphps(px,py);
*/

proc(0) = graphps(px,py);
    _protate = 1;
    _paxht = 0.25;
    _pnumht = 0.25;
    _ptitlht = 0.25;
    begwind;
    makewind(9,4,0,2.5,1);
    xy(px,py);
    endwind;
endp;