<*+M2EXTENSIONS *>
<*-CHECKDIV *>
<*-CHECKRANGE *>
<*-COVERFLOW *>
<*-IOVERFLOW*>
<*+STORAGE *>
<*-GENFRAME*>
<*+NOPTRALIAS*>
<*CPU="PENTIUM"*>
<* IF __GEN_C__ THEN *>
  <*-GENCTYPES*>
  <*+COMMENT*>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-GENDATE*>
  <*-LINENO*>
  <*-CHECKINDEX*>
  <*+GENCDIV*>
  <*-GENKRC*>
  <*+NOOPTIMIZE*>
  <*-GENSIZE*>
  <*-ASSERT*>
  <*-CHECKNIL*>
  <*-COVERFLOW*>
  <*-IOVERFLOW*>
  <*-CHECKRANGE*>
  <*-CHECKSET*>
  <*-CHECKDIV*>
  <*-GENCONSTENUM*>
<* ELSE *>
  <*-GENHISTORY*>
  <*+GENDEBUG*>
  <*+LINENO*>
  <*-CHECKINDEX*>
<* END *>

IMPLEMENTATION MODULE libsrtm; (* get altitude out of srtm files directory tree by oe5dxl *)

FROM osi IMPORT WrStr, WrStrLn, WrInt, FdValid, InvalidFd, File, Exists, Erase,
                RdBin, OpenWrite, WrBin, Close, OpenRead, DIRSEP, Seek,
                sqrt, sin, arctan, tan, ln, cos, pi, exp, WrFixed, floor, power,
                ALLOCATE, DEALLOCATE, time;
FROM SYSTEM IMPORT ADDRESS, FILL, ADR, INT16, INT8, CAST, MOVE, CARD8, CARD16, SHIFT, TSIZE;
FROM aprspos IMPORT posvalid, distance, azimuth, EARTH, RAD;
FROM aprsstr IMPORT postoloc, loctopos, POSITION, IntToStr, FixToStr, Append, Length,
                    posinval, TimeToStr, StrToFix, TIME, cleanfilename, InStr, Assign;


CONST 
      SRTMXY=3600;
      STRIPS=3;
      AGELEVELS=5;
      ATTRSUB=10000;             (* segment altitude to add metadata *)
      ATTRNEG=1000;              (* max under see level *)
      NOALT=32767;

TYPE
     SET8=SET OF [0..7];
     SET16=SET OF [0..15];
     FN=ARRAY[0..1023] OF CHAR;


    SRTMSTRIP=ARRAY[0..SRTMXY DIV STRIPS-1] OF INT16;
    pSRTMSTRIP=POINTER TO SRTMSTRIP;
    pSRTMTILE=POINTER TO SRTMTILE;
    SRTMTILE=RECORD
      typ   :CARD8;
      fd    :File;
      used  :ARRAY[0..STRIPS-1] OF ARRAY[0..SRTMXY-1] OF CARD8;
      strips:ARRAY[0..STRIPS-1] OF ARRAY[0..SRTMXY-1] OF pSRTMSTRIP;
    END;
    SRTMLAT=ARRAY[0..179] OF pSRTMTILE;
    pSRTMLAT=POINTER TO SRTMLAT;
    SRTMLONG=ARRAY[0..359] OF pSRTMLAT;
    SRTM30FD=ARRAY[0..8] OF ARRAY[0..3] OF RECORD
      fd:File;
      havefile:BOOLEAN;                                           (* have tried to open file *)
    END;
   

VAR
    srtmcache    :SRTMLONG;
    srtmmiss     :pSRTMTILE;                                (* cache no file info with pointer to here *) 
    srtm30fd     :SRTM30FD;                                 (* open srtm30 files *)
    errflag      :BOOLEAN;

(* === srtm lib *)
PROCEDURE opensrtm(t:CARD8; tlat, tlong:CARDINAL):File;
CONST SRTM3DIR="srtm3";
      SRTM1DIR="srtm1";
      SRTM30DIR="srtm30";

VAR s:ARRAY[0..20] OF CHAR;
    n, xi, yi, xd:CARDINAL;
    f:File;
    path:ARRAY[0..1024] OF CHAR;
BEGIN
  Assign(path, srtmdir);
  IF t=3 THEN Append(path, DIRSEP+SRTM3DIR+DIRSEP);
  ELSIF t=1 THEN Append(path, DIRSEP+SRTM1DIR+DIRSEP); END;
  IF t<=3 THEN
    IF tlat<90 THEN s[0]:="S"; n:=90-tlat ELSE s[0]:="N"; n:=tlat-90 END;
    s[1]:=CHR(n DIV 10+ORD("0"));
    s[2]:=CHR(n MOD 10+ORD("0"));
    IF tlong<180 THEN s[3]:="W"; n:=180-tlong ELSE s[3]:="E"; n:=tlong-180 END;
    s[4]:=CHR(n DIV 100+ORD("0"));
    s[5]:=CHR(n DIV 10 MOD 10+ORD("0"));
    s[6]:=CHR(n MOD 10+ORD("0"));
    s[7]:=".";
    s[8]:="h";
    s[9]:="g";
    s[10]:="t";
    s[11]:=0C;
    Append(path, s);
    RETURN OpenRead(path);

  ELSE
    Append(path, DIRSEP+SRTM30DIR+DIRSEP);
    xi:=tlong DIV 40;
    xd:=xi*40;
    IF xd<180 THEN s[0]:="W"; n:=180-xd ELSE s[0]:="E"; n:=xd-180 END;
    s[1]:=CHR(n DIV 100+ORD("0"));
    s[2]:=CHR(n DIV 10 MOD 10+ORD("0"));
    s[3]:=CHR(n MOD 10+ORD("0"));

    IF tlat>=130 THEN s[4]:="N"; s[5]:="9"; yi:=3;
    ELSIF tlat>=80 THEN s[4]:="N"; s[5]:="4"; yi:=2;
    ELSIF tlat>=30 THEN s[4]:="S"; s[5]:="1"; yi:=1;
    ELSE                                                   (* does not exist *)
      yi:=0; 
      srtm30fd[xi, yi].fd:=InvalidFd;
      srtm30fd[xi, yi].havefile:=TRUE;
    END;
--WrInt(xi, 4);WrInt(yi, 4); WrInt(ORD(srtm30fd[xi, yi].havefile), 4);
--WrInt(srtm30fd[xi, yi].fd, 5); WrStrLn(" x y have fd");

    IF srtm30fd[xi, yi].havefile THEN RETURN srtm30fd[xi, yi].fd END;

    s[6]:="0";
    s[7]:=".";
    s[8]:="D";
    s[9]:="E";
    s[10]:="M";
    s[11]:=0C;
    Append(path, s);
    f:=OpenRead(path);
    srtm30fd[xi, yi].fd:=f;
    srtm30fd[xi, yi].havefile:=TRUE;
    RETURN f
  END;
END opensrtm;


PROCEDURE purgesrtm(all:BOOLEAN);
CONST PURGEALL=10;
VAR xd, yd, x, y, asize, destmem:CARDINAL;
    pt:pSRTMTILE;
    pl:pSRTMLAT;
    pb:pSRTMSTRIP;
    empty:BOOLEAN; 

BEGIN
    destmem:=srtmmaxmem - srtmmaxmem DIV 10;
    FOR xd:=0 TO HIGH(srtmcache) DO
      pl:=srtmcache[xd]; 
      IF pl<>NIL THEN
        FOR yd:=0 TO HIGH(pl^) DO
          empty:=TRUE;
          pt:=pl^[yd];
          IF (pt<>NIL) & (pt<>srtmmiss) THEN
            WITH pt^ DO 
              FOR y:=0 TO SRTMXY-1 DO
                FOR x:=0 TO STRIPS-1 DO
                  pb:=strips[x][y]; 
                  IF pb<>NIL THEN
                    IF NOT all & (used[x][y]>0) THEN
                      DEC(used[x][y]);                                           (* increase age level *)
                      empty:=FALSE;
                    ELSE
                      asize:=SIZE(pb^);
                      IF typ>3 THEN asize:=SIZE(pb^) DIV 10 END;
                      DEALLOCATE(pb, asize);
                      DEC(srtmmem, asize);    
                      strips[x][y]:=NIL;
                      IF NOT all & NOT empty & (srtmmem<=destmem) THEN RETURN END; (* purged enough now *)

                    END;
                  END;
                END;
              END;
              IF all OR empty THEN
                IF (typ<=3) & (fd<>InvalidFd) THEN Close(fd) END;
                DEALLOCATE(pt, SIZE(pt^));
                DEC(srtmmem, SIZE(pt^));
                pl^[yd]:=NIL;
              END;
            END;
          END;
        END;
        IF all THEN
          DEALLOCATE(pl, SIZE(pl^));
          DEC(srtmmem, SIZE(pl^));
          srtmcache[xd]:=NIL;
        END;
      END;
    END;
    IF all THEN
      FOR x:=0 TO HIGH(srtm30fd) DO
        FOR y:=0 TO HIGH(srtm30fd[0]) DO
          IF (srtm30fd[x][y].havefile) & (srtm30fd[x][y].fd<>InvalidFd) THEN
            Close(srtm30fd[x][y].fd);
            srtm30fd[x][y].fd:=InvalidFd;
          END;
        END;
      END;
    END;
END purgesrtm;


PROCEDURE getsrtm1(ilat, ilong:CARDINAL; VAR div:CARDINAL):INTEGER;   (* 1 pixel altitude *)
VAR i, x, y, xx, ydeg, xdeg, rdsize:CARDINAL;
    seek:INTEGER;
    f:File;
    pt:pSRTMTILE;
    pb:pSRTMSTRIP;
    a:INTEGER; 
    t:CARD8;
BEGIN
  div:=1;
  ydeg:=ilat  DIV SRTMXY;
  xdeg:=ilong DIV SRTMXY;

  IF (xdeg>HIGH(srtmcache)) OR (ydeg>HIGH(srtmcache[0]^)) THEN RETURN NOALT END;
 
  IF srtmcache[xdeg]=NIL THEN                                      (* empty lat array *)
    ALLOCATE(srtmcache[xdeg], SIZE(srtmcache[0]^));;
    IF srtmcache[xdeg]=NIL THEN RETURN NOALT END;                  (* out of memory *)

    INC(srtmmem, SIZE(srtmcache[0]^));
    FILL(srtmcache[xdeg], 0C, SIZE(srtmcache[0]^));
  ELSIF srtmcache[xdeg]^[ydeg]=srtmmiss THEN RETURN NOALT END;     (* tile file not avaliable *)

  pt:=srtmcache[xdeg]^[ydeg];
  IF pt=NIL THEN
    t:=1;
    f:=opensrtm(t, ydeg, xdeg);
    IF f=InvalidFd THEN
      t:=3;
      f:=opensrtm(t, ydeg, xdeg);
      IF f=InvalidFd THEN
        t:=30;
        f:=opensrtm(t, ydeg, xdeg);
        IF f=InvalidFd THEN
          srtmcache[xdeg]^[ydeg]:=srtmmiss;
          RETURN NOALT

        END;

      END;
    END;
--INC(open);
    ALLOCATE(pt, SIZE(pt^));                                       (* a new 1x1 deg buffer *) 
    IF pt=NIL THEN RETURN NOALT END;

    INC(srtmmem, SIZE(pt^));
    WITH pt^ DO
      FILL(ADR(strips), 0C, SIZE(strips));
      FILL(ADR(used), 0C, SIZE(used));
      typ:=t;
      fd:=f;
    END;
    srtmcache[xdeg]^[ydeg]:=pt;
  END;

  WITH pt^ DO
    div:=typ; 
    IF typ=1 THEN
      y:=ilat  MOD SRTMXY;
      x:=ilong MOD SRTMXY; 
      xx:=x DIV (SRTMXY DIV STRIPS);
      x :=x MOD (SRTMXY DIV STRIPS);       
    ELSIF typ=3 THEN
      y:=(ilat  MOD SRTMXY) DIV 3;
      x:=(ilong MOD SRTMXY) DIV 3;
      xx:=0;
    ELSE
      y:=(ilat  MOD SRTMXY) DIV 30;
      x:=(ilong MOD SRTMXY) DIV 30;
      xx:=0;
    END;

    pb:=strips[xx][y];
    IF pb=NIL THEN
      IF typ=1 THEN 
        seek:=((SRTMXY-1)-y)*((SRTMXY+1)*2) + xx*((SRTMXY*2) DIV STRIPS);
        rdsize:=SIZE(pb^);
      ELSIF typ=3 THEN
        seek:=((SRTMXY DIV 3-1)-y)*((SRTMXY DIV 3+1)*2);
        rdsize:=SIZE(pb^);
      ELSE
        seek:=VAL(INTEGER, (xdeg MOD 40)*240) - VAL(INTEGER, (ilat DIV 30)*9600);
        IF ydeg>=130 THEN   INC(seek, 9600*((180*3600) DIV 30-1))
        ELSIF ydeg>=80 THEN INC(seek, 9600*((130*3600) DIV 30-1))
        ELSIF ydeg>=30 THEN INC(seek, 9600*((80 *3600) DIV 30-1))
        ELSE RETURN NOALT END;

        rdsize:=SIZE(pb^) DIV 10;                                      (* fill 1/10 buffer *) 
      END;
--INC(miss);

      ALLOCATE(pb, rdsize);
      IF pb=NIL THEN RETURN NOALT END;

      INC(srtmmem, rdsize);
      strips[xx][y]:=pb;

      Seek(fd, seek);
      IF RdBin(fd, pb^, rdsize)<>VAL(INTEGER, rdsize) THEN
        FOR i:=0 TO rdsize DIV 2-1 DO pb^[i]:=VAL(INT16, NOALT) END;
      ELSE 
(*
c-translator frissts nicht
        FOR i:=0 TO rdsize DIV 2-1 DO pb^[i]:=pb^[i]<<8 + pb^[i]>>8 END;  (* motorola format *)  
*)
        FOR i:=0 TO rdsize DIV 2-1 DO 
          pb^[i]:=CAST(INT16, SHIFT(CAST(SET16, pb^[i]), 8) + SHIFT(CAST(SET16, pb^[i]), -8));
        END;  (* motorola format *)
      END;


--ELSE INC(hit);
    END;
--    IF used[xx][y]<AGELEVELS THEN INC(used[xx][y]) END;

    used[xx][y]:=AGELEVELS;

    a:=pb^[x];
    IF (a>32000) OR (a<-32000) THEN RETURN NOALT END;

    RETURN a;   
  END;
END getsrtm1;

PROCEDURE qint(a,b,c,v:REAL):REAL;     (* spline interpolator *)
BEGIN
  a:=(a-b)*0.5;
  c:=(c-b)*0.5;

  RETURN b + (v*v+0.25)*(a+c) + v*(c-a)
END qint;

PROCEDURE qintd(a,b,c,v:REAL):REAL;     (* spline interpolator and normalvector *)
BEGIN RETURN v*(a+c-b*2.0) + (c-a)*0.5 END qintd;

(*
PROCEDURE qintd(a,b,c,v:REAL; VAR nv:REAL):REAL;     (* spline interpolator and normalvector *)
VAR ca, ac:REAL;
BEGIN
  a:=(a-b)*0.5;
  c:=(c-b)*0.5;
  ca:=c-a;
  ac:=a+c;
  nv:=2.0*v*ac + ca;
  RETURN b + (v*v+0.25)*ac + v*ca
END qintd;
*)

PROCEDURE chkmeta(VAR a:INTEGER):CARDINAL;
BEGIN RETURN VAL(CARDINAL, a+(4*ATTRSUB+ATTRNEG)) DIV ATTRSUB END chkmeta;

PROCEDURE resmeta(a:INTEGER):REAL;
BEGIN RETURN VAL(REAL, VAL(INTEGER, VAL(CARDINAL, a+(4*ATTRSUB+ATTRNEG)) MOD ATTRSUB) - ATTRNEG) END resmeta;

PROCEDURE am(ilat, ilong:CARDINAL):REAL;
VAR d:CARDINAL;
    a:INTEGER;
BEGIN
  a:=getsrtm1(ilat, ilong, d);
  IF a>=NOALT THEN errflag:=TRUE END;
  RETURN resmeta(a);
END am;

PROCEDURE int4(a0,a1,a2,a3, vx,vy:REAL):REAL;
BEGIN RETURN (a1*vx + a0*(1.0-vx))*(1.0-vy) + (a3*vx + a2*(1.0-vx))*vy END int4;
--BEGIN RETURN (a0*(1.0-vx) + a1*vx)*(1.0-vy) + (a2*(1.0-vx) + a3*vx)*vy END int4;

VAR a00, a0, a1, a2, a3, a4, a5, a6, a7, a8:REAL;
    at0,at1,at2,at3:CARD8;
    lastilat, lastilong, lastilat1, lastilong1, lastilat2, lastilong2, div:CARDINAL;
    lastatt0:CARD8;
 
PROCEDURE getsrtmlong(lat, long:LONGREAL; quality:CARDINAL; bicubic:BOOLEAN;
                      VAR resolution:REAL; VAR att0:CARD8; pmeta:pMETAINFO):REAL;
VAR ilat, ilong, div2, d:CARDINAL;
    ilatd, ilatdd, ilongd, ilongdd:CARDINAL;
--    a0, a1, a2, a3, a4, a5, a6, a7, a8:REAL;
    i0, i1, i2, i3:INTEGER;
    an1, vx, vy, vx5, vy5:REAL;
    dx0, dx1, dx2:REAL;
BEGIN
  INC(accesses);
  IF (srtmmaxmem>0) & (srtmmem>srtmmaxmem) THEN
    purgesrtm(FALSE);
  END;

  lat := (90.0*SRTMXY-0.5) + lat  * (SRTMXY/RAD);
  long:=(180.0*SRTMXY+0.5) + long * (SRTMXY/RAD);
  IF lat >=0.0 THEN ilat :=TRUNC(lat)  ELSE ilat:=0 END;
  IF long>=0.0 THEN ilong:=TRUNC(long) ELSE ilong:=0 END;

  IF (ilat<>lastilat) OR (ilong<>lastilong) THEN
    i0:=getsrtm1(ilat, ilong, div);
    lastatt0:=chkmeta(i0);
(*
  IF pmeta<>NIL THEN
    pmeta^.slantx:=0.0;
    pmeta^.slanty:=0.0;
    FILL(ADR(pmeta^.attrweights[0]), 0C, SIZE(pmeta^.attrweights));
  END;
*)
    IF i0>=NOALT THEN RETURN VAL(REAL,NOALT) END;
 
    a00:=resmeta(i0);
    lastilat:=ilat;
    lastilong:=ilong;
  ELSE INC(accessescached) END;
  att0:=lastatt0;

  IF pmeta<>NIL THEN quality:=0 END;              (* we always need 4 corners *)
  IF div=1 THEN
    resolution:=30.0;
    IF quality>29 THEN RETURN a00 END;

    vx:=long-LFLOAT(ilong);
    vy:=lat -LFLOAT(ilat);
  ELSIF div=3 THEN
    resolution:=90.0;
    IF quality>60 THEN RETURN a00 END;

    vx:=(long-LFLOAT((ilong DIV 3)*3))*(1.0/3.0);
    vy:=(lat -LFLOAT((ilat  DIV 3)*3))*(1.0/3.0);
  ELSE
    resolution:=900.0;
    IF quality>300 THEN RETURN a00 END;

    vx:=(long-LFLOAT((ilong DIV 30)*30))*(1.0/30.0);
    vy:=(lat -LFLOAT((ilat  DIV 30)*30))*(1.0/30.0);
  END;

  IF vx<0.5 THEN IF ilong>=div THEN DEC(ilong, div) END;
  ELSE vx:=vx-1.0 END;
  IF vy<0.5 THEN DEC(ilat,  div) ELSE vy:=vy-1.0 END;
  vx5:=vx+0.5;
  vy5:=vy+0.5;

  ilatd:=ilat+div;
  ilongd:=ilong+div;

  IF (ilat<>lastilat1) OR (ilong<>lastilong1) THEN
    i0:=getsrtm1(ilat, ilong, d);
    i1:=getsrtm1(ilat, ilongd, d);
    i2:=getsrtm1(ilatd, ilong, d);
    i3:=getsrtm1(ilatd, ilongd, d);
    IF (i0>=NOALT) OR (i1>=NOALT) OR (i2>=NOALT) OR (i3>=NOALT) THEN RETURN a00 END;

    a0:=resmeta(i0);
    a1:=resmeta(i1);
    a2:=resmeta(i2);
    a3:=resmeta(i3);
    at0:=chkmeta(i0);
    at1:=chkmeta(i1);
    at2:=chkmeta(i2);
    at3:=chkmeta(i3);
    lastilat1:=ilat;
    lastilong1:=ilong;
  END;
  
  IF (pmeta<>NIL) & pmeta^.aliasattr THEN         (* want attibute interpolation *)
    WITH pmeta^ DO
      FILL(ADR(attrweights[0]), 0C, SIZE(attrweights));
      attrweights[at0]:=attrweights[at0]+(1.0-vx5)*(1.0-vy5);
      attrweights[at1]:=attrweights[at1]+vx5*(1.0-vy5);
      attrweights[at2]:=attrweights[at2]+(1.0-vx5)*vy5;
      attrweights[at3]:=attrweights[at3]+vx5*vy5;
    END;
  END;

--interpolate 4 dots
  IF NOT bicubic THEN                                        (* bilinear *)

    IF (pmeta<>NIL) & pmeta^.withslant THEN                  (* bilinear slants *)
      pmeta^.slantx:=((a1 - a0)*(1.0-vy5) + (a3 - a2)*vy5)/resolution;
      pmeta^.slanty:=((a2 - a0)*(1.0-vx5) + (a3 - a1)*vx5)/resolution;
    END;

    RETURN int4(a0,a1,a2,a3,vx5,vy5);                        (* wighted median of 4 pixels*)
    
  ELSE                                                       (* biqubic *)
--    IF vx<0.5 THEN IF ilong>=div THEN DEC(ilong, div) END;
--    ELSE vx:=vx-1.0 END;
--    IF vy<0.5 THEN DEC(ilat,  div) ELSE vy:=vy-1.0 END;
    div2:=div*2;
    ilatdd:=ilat+div2;
    ilongdd:=ilong+div2;
    errflag:=FALSE;
    IF (pmeta<>NIL) & pmeta^.withslant THEN                  (* cubic interpolate and slants *)
      IF (ilat<>lastilat2) OR (ilong<>lastilong2) THEN
        a4:=am(ilat,   ilongdd);
        a5:=am(ilatd,  ilongdd);
        a6:=am(ilatdd, ilong);
        a7:=am(ilatdd, ilongd);
        a8:=am(ilatdd, ilongdd);
        lastilat2:=ilat;
        lastilong2:=ilong;
      END; 
      dx0:=qint(a0, a1, a4, vx);
      dx1:=qint(a2, a3, a5, vx);
      dx2:=qint(a6, a7, a8, vx);
      pmeta^.slanty:=qintd(dx0, dx1, dx2, vy)/resolution;

      dx0:=qint(a0, a2, a6, vy);
      dx1:=qint(a1, a3, a7, vy);
      dx2:=qint(a4, a5, a8, vy);
      pmeta^.slantx:=qintd(dx0, dx1, dx2, vx)/resolution;
      an1:=qint(dx0, dx1, dx2, vx);
      IF errflag THEN RETURN a00 ELSE RETURN an1 END; 

    ELSE an1:=qint(
      qint(a0, a1, am(ilat,   ilongdd), vx),
      qint(a2, a3, am(ilatd,  ilongdd), vx),
      qint(am(ilatdd, ilong), am(ilatdd, ilongd), am(ilatdd, ilongdd), vx), vy);
      IF errflag THEN RETURN a00 ELSE RETURN an1 END;

    END;

  END;

END getsrtmlong;

(*
PROCEDURE getsrtm(pos:POSITION; quality:CARDINAL; VAR resolution:REAL):REAL;
VAR attr:CARD8;
BEGIN RETURN getsrtmlong(pos.lat, pos.long, quality, FALSE, resolution, attr, NIL) END getsrtm;
*)
(*
PROCEDURE getsrtm(pos:POSITION; quality:CARDINAL; VAR resolution:REAL):REAL;
VAR attr:CARD8;
    a:REAL;
BEGIN
  a:=getsrtmlong(pos.lat, pos.long, quality, bicubic, resolution, attr, NIL);
  IF attr=3 THEN a:=32000.0 END;
  RETURN a
END getsrtm;
*)

PROCEDURE getsrtm(pos:POSITION; quality:CARDINAL; VAR resolution:REAL):REAL;
VAR attr:CARD8;
    a:REAL;
    at:METAINFO;
BEGIN
  at.withslant:=TRUE;
  at.aliasattr:=TRUE;

  a:=getsrtmlong(pos.lat, pos.long, quality, FALSE, resolution, attr, ADR(at));
--  IF (attr=3) OR (at.attr2=4) & ((attr=3) = (at.attrweight<0.5)) THEN a:=32000.0 END;
  
--IF attr=2 THEN a:=32000.0 END;
--  IF at.attrweights[2]>0.5 THEN a:=32000.0 END;
--  a:=1000.0+at.slanty*2.0;
  RETURN a
END getsrtm;


PROCEDURE initsrtm;
VAR xi, yi:CARDINAL;
BEGIN

--WrStrLn(" inits");

--open:=0; miss:=0; hit:=0;
  FILL(ADR(srtmcache), 0C, SIZE(srtmcache));
  FOR xi:=0 TO HIGH(srtm30fd) DO
    FOR yi:=0 TO HIGH(srtm30fd[0]) DO srtm30fd[xi][yi].havefile:=FALSE END;
  END;
END initsrtm;

PROCEDURE closesrtmfile;
BEGIN

--WrStrLn(" closes");

  purgesrtm(TRUE);
--WrInt(open, 10); WrInt(miss, 10); WrInt(hit, 10); WrStrLn(" open miss hit");
  initsrtm;
END closesrtmfile;


(* === srtm lib *)


(* === geoid lib *)

CONST GEOIDFN="WW15MGH.DAC";
      RESOL=4;                      (* 1/deg *)
      LONGS=360*RESOL;              (* values in file around a latitude *) 

PROCEDURE rdgeoid(fd:File; lat, long:INTEGER; VAR ok:BOOLEAN):REAL;
VAR b:ARRAY[0..1] OF CHAR;
    n:INTEGER;
BEGIN
  Seek(fd, (lat*LONGS + long)*2);
  IF RdBin(fd, b, 2)<>2 THEN                      (* no data in file *)      
    n:=0;
    ok:=FALSE;
  ELSE
    n:=ORD(b[1])+ORD(b[0])*256;
    IF n>=32768 THEN DEC(n,65536) END;
  END;
  RETURN VAL(REAL, n)
END rdgeoid;


PROCEDURE egm96(pos:POSITION; VAR ok:BOOLEAN):REAL;   (* read and interpolate geoid correction from datafile *)
VAR ilat, ilong    :INTEGER;
    fd             :File;
    flat, flong, g :REAL;
    path           :ARRAY[0..1024] OF CHAR;

BEGIN
  
  Assign(path, srtmdir);
  IF path<>"" THEN Append(path, DIRSEP) END;
  Append(path, GEOIDFN);
  fd:=OpenRead(path);
  ok:=fd<>InvalidFd;

  g:=0.0;
  IF ok THEN
    pos.lat:=90.0-pos.lat*(1.0/RAD);
    pos.long:=pos.long*(1.0/RAD);
    IF pos.long<0.0 THEN pos.long:=360.0 + pos.long END;
    pos.lat:=pos.lat*RESOL;
    pos.long:=pos.long*RESOL;
    IF pos.lat >=0.0 THEN ilat :=TRUNC(pos.lat)  ELSE ilat :=0 END;
    IF pos.long>=0.0 THEN ilong:=TRUNC(pos.long) ELSE ilong:=0 END;
    flat:=pos.lat-FLOAT(ilat);
    flong:=pos.long-FLOAT(ilong);
    g:=(( rdgeoid(fd, ilat, ilong  , ok)*(1.0-flat)
        + rdgeoid(fd, ilat+1, ilong  , ok)*flat)*(1.0-flong)
        +(rdgeoid(fd, ilat, ilong+1, ok)*(1.0-flat)
        + rdgeoid(fd, ilat+1, ilong+1, ok)*flat)*flong)*0.01;
    Close(fd);
  END;
  RETURN g
END egm96;


BEGIN
  ALLOCATE(srtmmiss, 1);                        (* make empty tile for fast nofile hint *)
  initsrtm;
  srtmmem:=0;
  srtmmaxmem:=0;                                (* 0 no auto purge *)
  srtmdir:="";
  bicubic:=FALSE;
  lastilat:=MAX(CARDINAL);
  lastilong:=MAX(CARDINAL);
  lastilat1:=MAX(CARDINAL);
  lastilong1:=MAX(CARDINAL);
  lastilat2:=MAX(CARDINAL);
  lastilong2:=MAX(CARDINAL);
  accessescached:=0;
  accesses:=0;
END libsrtm.