(*
gcc  -o panorama Lib.o aprspos.o aprsstr.o flush.o libsrtm.o osi.o pngwrite.o panorama.o /usr/local/xds/lib/x86/libts.a /usr/local/xds/lib/x86/libxds.a -lm libz.a filesize.o imagetext.o pngread.o libpng.a libz.a libjpeg.a
*)
<*+M2EXTENSIONS *>
<*-CHECKDIV *>
<*-CHECKRANGE *>
<*-GENFRAME*>
<*-COVERFLOW *>
<*-IOVERFLOW*>
<*-NOPTRALIAS*>
<*-DOREORDER*>
<*-PROCINLINE*>
<*-GENPTRINIT*>
<*+STORAGE *>
<*CPU="PENTIUM"*>
<* IF __GEN_C__ THEN *>
  <*+M2EXTENSIONS *>
  <*+STORAGE *>
  <*-GENCTYPES*>
  <*+COMMENT*>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-GENDATE*>
  <*-LINENO*>
  <*-CHECKINDEX*>
  <*-CHECKDINDEX*>
  <*+GENCDIV*>
  <*-GENKRC*>
  <*+NOOPTIMIZE*>
  <*-GENSIZE*>
  <*-ASSERT*>
  <*-CHECKNIL*>
  <*-COVERFLOW*>
  <*-IOVERFLOW*>
  <*-CHECKRANGE*>
  <*-CHECKSET*>
  <*-CHECKDIV*>
  <*-GENCONSTENUM*>
  <*-ASSERT*>
<* ELSE *>
  <*+GENHISTORY*>
  <*+GENDEBUG*>
  <*+LINENO*>
  <*+CHECKNIL*>
  <*+CHECKINDEX*>
  <*+CHECKDINDEX*>
  <*HEAPLIMIT="1000000000"*>

<* END *>

MODULE panorama; (* make panorama image with srtm data by oe5dxl *)

FROM SYSTEM IMPORT FILL, ADR, INT16, INT8, CAST, CARD8, CARD16, SHIFT;

FROM osi IMPORT Werr, WrStr, WrStrLn, WrInt, File, OpenWrite, Close, WrBin, 
     WrFixed, ALLOCATE, DEALLOCATE, NextArg, RdBin, FdValid, OpenRead, Flush;

FROM libsrtm IMPORT getsrtm, srtmdir, closesrtmfile, srtmmaxmem, getsrtmlong,
                    METAINFO, pMETAINFO, ATTRWATER, ATTRWOOD, ATTRURBAN,
                    accessescached, accesses;

FROM aprsstr IMPORT POSITION, IntToStr, Append, Assign, DateToStr, FixToStr, Length,
     loctopos,  StrToFix, StrToInt, StrToCard, TIME, InStr;

IMPORT pngwrite, pngread, osi;
FROM imagetext IMPORT pIMAGE, writestr, strsize, fontsizex, fontsizey;
FROM math IMPORT sqrt, sin, cos, tan, atan, exp, log, pow, floor, fmod, asin;

CONST
  TESTX=2;
  TESTY=359;
      PI=3.1415926535;
      PI2=PI*2.0;
      RAD=PI2/360.0;
      EARTH=6370.0;
      NOALT=20000;
      NL=12C;
      LF=12C;   
      REFRACC=0.0000000785;
      MAXSYMSIZE=48;                                        (* max input symbol size *)
      SUNSIZE=0.55*0.5*RAD;
      SHOWTREEDIST=20;                                      (* show trees m from cam *)
      ZEROTREEALT=2600;
      FULLTREEALT=1200;
      TREERASTERSIZE=8.0;                                   (* default m full treesize *)
      STAMTHICKNESS=0.015;
      TREEFORM=0.9; (* 0.9 *)
      STAMLEN=0.35;
      MINTREESIZE=3;
      NOISESIZE=1024;
      MAXHIST=250;
      MINHIST=10;
      MAXLUM=16383;
      TREERAND=3;                                           (* tree position noise *)
      LUMMUL=8.0;                                           (* real to cardinal colour values in main image *)
      ALTICOLORANG=1000;                                    (* colourized altitude *)
      WATERMIRRORGAP=0.2;                                   (* m over water start ray path *)
      WATERSURFACEDUST=0.6;
      WAVLEN=5000000.0;

      TREENOTCH=0.3;                                       (* random tree kill *)

      cFOGLUMr=900.0;          (* low elevation heaven light *)
      cFOGLUMg=1250.0;
      cFOGLUMb=2300.0;

      SUNLUM=5000.0;
      cMAXCHROMr=1.0;          (* sun passed thu air *)
      cMAXCHROMg=0.9;
      cMAXCHROMb=0.6;          (*0.8*)

      cFARCHROMr=1.0;          (* sun in space *)
      cFARCHROMg=1.0;
      cFARCHROMb=1.0;

      cMEADOWr=310.0*0.8;
      cMEADOWg=430.0*0.8;
      cMEADOWb=215.0*0.8;      (* 235 *)
      cMEADOWVOIDr=458.0*0.9;
      cMEADOWVOIDg=458.0*0.9;
      cMEADOWVOIDb=429.0*0.9;
      cCOASTSANDr=380.0;
      cCOASTSANDg=450.0;
      cCOASTSANDb=450.0;

      cSNOWr=840.0;            (*940 930 920 *)
      cSNOWg=830.0;
      cSNOWb=820.0;
      cSNOWSHADOWr=79.0;       (*100 145 280*)
      cSNOWSHADOWg=115.0;
      cSNOWSHADOWb=220.0;

      cSTAMr=140.0*0.7;            (* 292 271 249 *)
      cSTAMg=130.0*0.7;
      cSTAMb=120.0*0.7;
      cSTAMVAR=15.0;
      cTREEr=210.0*0.72;            (* 200 245 190  *0.85 *)
      cTREEg=270.0*0.72;
      cTREEb=225.0*0.72;

      cCITYCOL=400.0;
      cCITYMEADOWr=258.0;
      cCITYMEADOWg=366.0;
      cCITYMEADOWb=292.0;
      cCITYSTREETr=249.0;
      cCITYSTREETg=282.0;
      cCITYSTREETb=282.0;

      cCITYDUSTr=0.18;
      cCITYDUSTg=0.15;
      cCITYDUSTb=0.1;

      cROCKr=418;
      cROCKg=458;
      cROCKb=448;

      cUNDERWATERr=360.0*0.5;     (* 360 360 330 * 0.6 *)
      cUNDERWATERg=320.0*0.5;
      cUNDERWATERb=300.0*0.5;

      cWATERALIASr=120.0;     (* basic watercolour 33 81 99 *0.9 *)
      cWATERALIASg=170.0;
      cWATERALIASb=220.0;

      cHEAVENr=150.0;
      cHEAVENg=300.0;
      cHEAVENb=3500.0;

      cLAMPHr=0.3;
      cLAMPHg=0.24;
      cLAMPHb=0.2;
      cLAMPWr=0.3;
      cLAMPWg=0.3;
      cLAMPWb=0.3;
      cLAMPTr=0.4;
      cLAMPTg=0.1;
      cLAMPTb=0.0;

TYPE
     SET8=SET OF [0..7];
     SET16=SET OF [0..15];
     SET32=SET OF [0..31];
     COL=RECORD r,g,b:REAL END;
     CCOL=RECORD r,g,b:CARDINAL END;
     SCREENTYPE=(eCURVED, eFLAT, eCAMERA);
     POSITIONL=RECORD
                   long,
                   lat:LONGREAL;
                 END;

     PANOWIN=RECORD
       flatscreen:SCREENTYPE;
       eye, horizon:POSITIONL;
       eyealt, eyeog, angle, elevation, yzoom, refract:REAL;
       hx, hy, mx, my:INTEGER;                (* mouse pos on panwin and fullwin *)
       minlum, maxlum:CARDINAL;               (* hint for text contrast *)
       sunazi, sunele,
       hseg, dist:REAL;                  (* poi fast preselect circular segment *)
       image:pIMAGE;
     END;

     PIX8A=RECORD r8, g8, b8, alpha:CARD8 END;                  (* colour pixel with alpha *)
     SYMBOL=ARRAY[0..MAXSYMSIZE-1] OF ARRAY[0..MAXSYMSIZE-1] OF PIX8A;

     pPOI=POINTER TO POI;
     POI=RECORD
       next:pPOI;
       r,g,b,alpha:CARDINAL;
       symbol, woodsymbol:SYMBOL;
       xs,ys:CARDINAL;
       altIsOG, withwoodicon:BOOLEAN;
       prio:CARD8;
       poifn, iconfn, iconwoodfn, iconhint:ARRAY[0..1023] OF CHAR;
     END;

     pPOISTORE=POINTER TO POISTORE;
     POISTORE=RECORD
       next:pPOISTORE;
       pp:pPOI;
       x, y:REAL;
       prio:CARDINAL;
       dist:REAL;
       thrutree:BOOLEAN;
       textposx:INTEGER;
       text:ARRAY[0..50] OF CHAR;
     END; 

     pWAYPOINTS=POINTER TO WAYPOINTS;
     WAYPOINTS=RECORD
       wps:ARRAY[0..63] OF RECORD pos:POSITIONL; alt:REAL END;
       next:pWAYPOINTS;
     END; 


     pPOIMAP=POINTER TO POIMAP;                         (* store image columns with poi text *)
     POIMAP=ARRAY OF SET32;

     pTRACK=POINTER TO TRACK;
     TRACK=RECORD
       next:pTRACK;
       r,g,b, altmode:CARDINAL;
       width, fd:INTEGER;
       cnt:CARDINAL;
       fn:ARRAY[0..99] OF CHAR;
     END;

     NORMVEC=RECORD x,y,z:REAL END;
     LVEC=RECORD x,y,z:LONGREAL END;
     ROTVEC=RECORD sx,cx,sy,cy:LONGREAL END;
     RGB=RECORD r,g,b:REAL END;
     GLANCESET=RECORD sunvec, eyevec, surfacevec:NORMVEC END;

     TREEHINT=RECORD
                hith,
                hitx, hity,
                sx,sy,
                peakh,
                areapeak,
                regionpeak,
                hitradius,
                stamtoeye,
                og,
                passdist,
                aliasm,
                missspace:REAL;
                nx, ny:INTEGER;
                treepass,
                isstam:BOOLEAN;
                tslant,
                varh,
                rooth:REAL;
                meta:METAINFO;
                nvec:NORMVEC;
                water, wood, urban:CARD8;
     END;
     TREEHINTS=ARRAY[FALSE..TRUE] OF TREEHINT;

     NOISEMAP=ARRAY[0..NOISESIZE-1] OF ARRAY[0..NOISESIZE-1] OF REAL;
     ALTCOLTAB=RECORD
                 c:ARRAY[0..1023] OF CCOL;
                 tag:ARRAY[0..5] OF CCOL;
                 len:CARDINAL;
               END; 

     pDROPSHADOW=POINTER TO DROPSHADOW;
     DROPSHADOW=RECORD
                  next:pDROPSHADOW;
                  dist,
                  alt:LONGREAL;
                  colum:INTEGER;
                END;

     DDSTACK=RECORD
               startdist:REAL;
               stoppalt:REAL;
             END;

VAR
    image       :pIMAGE;
    poifiles    :pPOI;
    trackfiles  :pTRACK;
    dateh, timeh,      
    imagefn,
    csvfn, colfn:ARRAY[0..1023] OF CHAR;
    xsize, ysize:INTEGER; 
    posa, posb,
    teye        :POSITIONL;
    refraction, igamma, elev, angle, yzoom, wavelength, waterslant:REAL;
    gammatab    :ARRAY[0..MAXLUM] OF CHAR;
    millisec,
    altcolstep,
    rastercsv,
    poisize,
    scalesize,
    whitelim,
    blacklim,
    poilabel,
    fonttypfoot,
    fonttyp     :CARDINAL;
    watercol,
    heavencol   :COL;
    desertcol   :RGB;
    dropshadows,
    alticolour,
    defaultocean,
    verb,
    flatwater   :BOOLEAN;
    flatscreen  :SCREENTYPE;
    treeraster,
    coslat,
    treerasterlat,
    treerasterlong,
    treescramblat,
    treescramblong,
    maxmountain,
    rndseed:    LONGREAL;
    surfacedust,
    waterglance,
    watertransparency,
    treerastersize,
    abscolmin,
    abscolmax,
    alticolmin,
    alticolmax,
    zerotreealt,
    fulltreealt,
    treeshrink,
    camog,
    camalt,
    sunaz, sunel,
    moonazim,moonelev, moondist,
    snowalt,
    slowaltslant,
    newsnow,
    waveamp,
    dusthighpass, 
    treesize, 
    botanicslant,
    ignorblack,
    ignorwhite, 
    showstamdist,
    videolat,
    logfilm,
    urbanoff,
    sandgranularity,
    minpoialt   :REAL;
    csvfd       :INTEGER;
    csvvec      :RECORD
                   act:CARDINAL;
                   v:ARRAY[0..999] OF RECORD
                     p:POSITION;
                     a:REAL;
                   END;
                 END;
    poirect0, poirect1:POSITION;      (* poi fast preselect rectangle *)

    CRCL, CRCH : ARRAY[0..255] OF CARD8;
    CRC: ARRAY[0..65535] OF CARD16;
    noisemap:NOISEMAP;
    noisevec:ARRAY[0..8191] OF REAL;
    altcoltab:ALTCOLTAB;
    freeshadow,
    newshadow:pDROPSHADOW;
    poistore:pPOISTORE;
    poimap:pPOIMAP;

debugxi,debugyi:INTEGER;
debugx, debugy, debugz, debugg:REAL;

PROCEDURE StrToFixL(VAR x:LONGREAL; s-:ARRAY OF CHAR):BOOLEAN;
CONST M=MAX(LONGREAL)*0.01;
VAR i:CARDINAL;
    p:LONGREAL;
    neg, ok:BOOLEAN;
BEGIN
  i:=0;
  neg:=FALSE;
  ok:=FALSE;
  WHILE (i<=HIGH(s)) & (s[i]=" ") DO INC(i) END;
  IF s[i]="-" THEN neg:=TRUE; INC(i) END;
  x:=0.0;
  WHILE (i<=HIGH(s)) & (s[i]>="0") & (s[i]<="9") DO
    IF x<M THEN x:=x*10.0 + LFLOAT(ORD(s[i])-ORD("0")) END;
    INC(i);
    ok:=TRUE;
  END;
  IF (i<=HIGH(s)) & (s[i]=".") THEN
    INC(i);
    p:=0.1;
    WHILE (i<=HIGH(s)) & (s[i]>="0") & (s[i]<="9") DO
      x:=x + p*LFLOAT(ORD(s[i])-ORD("0"));
      p:=p*0.1;
      INC(i);
      ok:=TRUE;
    END;
  END;
  IF neg THEN x:=-x END;
  RETURN ok & ((i>HIGH(s)) OR (s[i]=0C))
END StrToFixL;



PROCEDURE Error(text:ARRAY OF CHAR);
BEGIN
  Werr(text); Werr(" error abort"+LF);
  HALT
END Error;

PROCEDURE hex(n:CARDINAL):CHAR;
BEGIN
  n:=n MOD 16;
  IF n<=9 THEN RETURN CHR(n+ORD("0")) END;
  RETURN CHR(n+(ORD("A")-10))
END hex;

PROCEDURE WrHex(x, digits, len:CARDINAL);
VAR i:CARDINAL;
    s:ARRAY[0..255] OF CHAR;
BEGIN
  IF digits>HIGH(s) THEN digits:=HIGH(s) END;
  i:=digits;
  WHILE (i<len) & (i<HIGH(s)) DO
    s[i]:=" ";
    INC(i);
  END;
  s[i]:=0C;
  WHILE digits>0 DO
    DEC(digits);
    s[digits]:=hex(x);
    x:=x DIV 16;
  END;
  WrStr(s); 
END WrHex;




PROCEDURE Random():LONGREAL;
BEGIN
  rndseed:=rndseed + PI;
  rndseed:=rndseed*rndseed;
  rndseed:=rndseed*rndseed;
  rndseed:=rndseed*rndseed;
  rndseed:=rndseed - LFLOAT(TRUNC(rndseed));
  RETURN rndseed
END Random;

PROCEDURE truncc(r:REAL):CARDINAL;
BEGIN
  IF r<=0.0 THEN RETURN 0 ELSIF r>=2000000000.0 THEN RETURN 2000000000 ELSE RETURN TRUNC(r) END;

<* IF NOT __GEN_C__ THEN *>

EXCEPT
  WrStrLn("error in TRUNC() ");
  RETURN 0
<* END *>
END truncc;

PROCEDURE Gencrctab;
CONST POLINOM=08408H;
VAR i,crc,c:CARDINAL;
BEGIN
  FOR c:=0 TO 255 DO
    crc:=255-c;
    FOR i:=0 TO 7 DO
      IF ODD(crc) THEN crc:=CAST(CARDINAL, CAST(BITSET, ASH(crc, -1))/CAST(BITSET,POLINOM))
      ELSE crc:=ASH(crc, -1) END;
    END;
    CRCL[c]:=crc;
    CRCH[c]:=255 - ASH(crc, -8);
  END;
  crc:=1;
  FOR i:=0 TO HIGH(CRC) DO 
    IF ODD(crc) THEN crc:=CAST(CARDINAL, CAST(BITSET, ASH(crc, -1))/CAST(BITSET,POLINOM))
    ELSE crc:=ASH(crc, -1) END;
    CRC[i]:=crc;
  END;
END Gencrctab;


PROCEDURE frac(r:LONGREAL):REAL;
BEGIN RETURN r-floor(r) END frac;

PROCEDURE fracint(r:LONGREAL; VAR i:INTEGER):REAL;
VAR ri:LONGREAL;
BEGIN
  ri:=floor(r);
  i:=VAL(INTEGER, ri);
  RETURN r-ri;
END fracint;

PROCEDURE min(a,b:REAL):REAL; BEGIN IF a<b THEN RETURN a END; RETURN b END min;
PROCEDURE max(a,b:REAL):REAL; BEGIN IF a>b THEN RETURN a END; RETURN b END max;
PROCEDURE imin(a,b:INTEGER):INTEGER; BEGIN IF a<b THEN RETURN a END; RETURN b END imin;
PROCEDURE imax(a,b:INTEGER):INTEGER; BEGIN IF a>b THEN RETURN a END; RETURN b END imax;
PROCEDURE sqr(x:REAL):REAL; BEGIN RETURN x*x END sqr;
PROCEDURE plus(x:REAL):REAL; BEGIN IF x>=0.0 THEN RETURN x END; RETURN 0.0 END plus;

PROCEDURE satan(x:REAL):REAL;   (* satturating tan at +-pi/2 *)
CONST LIM=PI*0.5-0.0001;
BEGIN
  IF x>=LIM THEN x:=LIM
  ELSIF x<=-LIM THEN x:=-LIM END;
  RETURN tan(x)
END satan;


PROCEDURE one(x:REAL):REAL;
BEGIN IF x<0.0 THEN x:=0.0 ELSIF x>1.0 THEN x:=1.0 END; RETURN x END one;

PROCEDURE atang2(x, y:REAL):REAL;
VAR w:REAL;
BEGIN
  IF ABS(x)>ABS(y) THEN
    w:=osi.arctan(y/x);
    IF x<0.0 THEN
      IF y>0.0 THEN w:=PI+w ELSE w:=w-PI END;
    END;
  ELSIF y<>0.0 THEN
    w:=(PI*0.5-osi.arctan(x/y));
    IF y<0.0 THEN w:=w-PI END;
  ELSE w:=0.0 END;
  RETURN w
END atang2;

PROCEDURE atang2l(x, y:LONGREAL):LONGREAL;
VAR w:LONGREAL;
BEGIN
  IF ABS(x)>ABS(y) THEN
    w:=atan(y/x);
    IF x<0.0 THEN
      IF y>0.0 THEN w:=PI+w ELSE w:=w-PI END;
    END;
  ELSIF y<>0.0 THEN
    w:=(PI*0.5-atan(x/y));
    IF y<0.0 THEN w:=w-PI END;
  ELSE w:=0.0 END;
  RETURN w
END atang2l;

PROCEDURE posvalidl(pos-:POSITIONL):BOOLEAN;
BEGIN RETURN (pos.lat<>0.0) OR (pos.long<>0.0) END posvalidl;

PROCEDURE distancel(home, dist:POSITIONL):REAL;
VAR x, y:REAL;
BEGIN
  x:=ABS(dist.long-home.long);
  y:=ABS(dist.lat -home.lat);
  IF (x=0.0) & (y=0.0) THEN RETURN 0.0
  ELSIF x+y<0.04 THEN                      (* near *)
    x:=x*cos((home.lat+dist.lat)*0.5);
    RETURN EARTH*sqrt(x*x+y*y)
  ELSE                                     (* far *)
    IF x>PI2 THEN x:=PI2-x END;
    x:=sin(dist.lat)*sin(home.lat) + cos(dist.lat)*cos(home.lat)*cos(x);
    IF ABS(x)>=1.0 THEN RETURN 0.0 END;
    RETURN EARTH*osi.arccos(x)
  END;
END distancel;


PROCEDURE azimuthl(home,dist:POSITIONL):REAL;
VAR h, ldiff:REAL;
BEGIN
  ldiff:=dist.long-home.long;
  IF (ldiff=0.0) OR (cos(home.lat)=0.0) OR (cos(dist.lat)=0.0) THEN
    IF home.lat<=dist.lat THEN RETURN 0.0 ELSE RETURN 180.0 END;
  ELSE
    h:=360.0/PI2*osi.arctan(cos(home.lat)*((tan(home.lat)*cos(ldiff)-tan(dist.lat))/sin(ldiff)));
    IF ldiff<0.0 THEN RETURN h+270.0 ELSE RETURN h+90.0 END;
  END;
END azimuthl;





PROCEDURE makenoise;
VAR x, y:CARDINAL;
    v:REAL;
PROCEDURE swap(w:CARDINAL):CARDINAL;
BEGIN RETURN w MOD 256*256+w DIV 256 END swap;

BEGIN
  FILL(ADR(noisemap), 0C, SIZE(noisemap));
  FOR y:=0 TO NOISESIZE DO
    FOR x:=0 TO NOISESIZE-1 DO
      v:=Random()-0.5;
      noisemap[y MOD NOISESIZE][x]:=noisemap[y MOD NOISESIZE][x]*0.7 + v;
      noisemap[y MOD NOISESIZE][x]:=noisemap[y MOD NOISESIZE][x]*0.3333 + noisemap[(y+(NOISESIZE-1)) MOD NOISESIZE][x]*0.3333 + noisemap[(y+1) MOD NOISESIZE][x]*0.3333;
    END;
  END;

  FILL(ADR(noisevec), 0C, SIZE(noisevec));
  FOR x:=0 TO HIGH(CRC) DO
    noisevec[x MOD (HIGH(noisevec)+1)]:=noisevec[x MOD (HIGH(noisevec)+1)]*0.7 + FLOAT(CRC[x])*(1.0/65536.0)-0.5;
    noisevec[x MOD (HIGH(noisevec)+1)]:=noisevec[(x+HIGH(noisevec)) MOD (HIGH(noisevec)+1)]*0.3333 + noisevec[x MOD (HIGH(noisevec)+1)]*0.3333 + noisevec[(x+1) MOD (HIGH(noisevec)+1)]*0.3333;
  END;

  v:=0.0;
  FOR x:=1 TO HIGH(noisevec) DO
    v:=noisevec[x]-noisevec[x-1];
    IF v>0.1 THEN v:=0.1 ELSIF v<-0.1 THEN v:=-0.1 END;
    noisevec[x]:=noisevec[x-1]+v;    
(*
    IF noisevec[x]>0.0 THEN noisevec[x]:=sqrt(noisevec[x])
    ELSIF noisevec[x]<0.0 THEN noisevec[x]:=-sqrt(-noisevec[x]) END;
    v:=v+(noisevec[x]-v)*0.2;
    noisevec[x]:=noisevec[x]-v;
*)
  END;
 

  v:=0.0;                                      (* normalize level *)
  FOR x:=0 TO HIGH(noisemap) DO
    FOR y:=0 TO HIGH(noisemap) DO
      v:=v+ABS(noisemap[y][x]);
    END;
  END;
  v:=sqr(FLOAT(HIGH(noisemap)+1))/v;
  FOR x:=0 TO HIGH(noisemap) DO
    FOR y:=0 TO HIGH(noisemap) DO
      noisemap[y][x]:=noisemap[y][x]*v;
    END;
  END;

--  WrFixed(v/sqr(FLOAT(HIGH(noisemap)+1)), 3, 10); WrStrLn("");

END makenoise;


(*
PROCEDURE rndx(x:LONGREAL):REAL;
BEGIN
  x:=x*10.0;
  RETURN (sin(x*072030.0)*0.0020-sin(x*1721560.0)*0.0014+sin(x*5303500.0)*0.0033
         -sin(x*787666.0)*0.0044+sin(x*1287666.0)*0.0051-sin(x*2355636.0)*0.0065)*5.5;
END rndx;
PROCEDURE rndy(x:LONGREAL):REAL;
BEGIN
  x:=x*10.0;
  RETURN (sin(x*104830.0)*0.0014-sin(x*2741560.0)*0.0026-sin(x*6373800.0)*0.0035
         +sin(x*986626.0)*0.0032-sin(x*2247676.0)*0.0053+sin(x*3254666.0)*0.0047)*5.5;
END rndy;
*)

PROCEDURE rndx(x:LONGREAL):REAL;
VAR a,b,i:CARDINAL;

BEGIN
  a:=VAL(INTEGER, x*1000000.0);
  b:=0;
  FOR i:=0 TO 4 DO
    INC(b, b+ORD(ODD(a)));
    a:=a DIV 2;
  END;
  x:=FLOAT(b);
  IF ODD(a) THEN x:=-x END;
  RETURN x*0.005
END rndx;

PROCEDURE rndy(x:LONGREAL):REAL;
VAR a,b,i:CARDINAL;
BEGIN
  a:=VAL(INTEGER, x*1000000.0);
  b:=0;
  FOR i:=0 TO 4 DO
    INC(b, b+ORD(ODD(a)));
    a:=a DIV 2;
  END;
  x:=FLOAT(b);
  IF ODD(a) THEN x:=-x END;
  RETURN x*0.006;
END rndy;





PROCEDURE posinval(VAR pos:POSITIONL);
BEGIN pos.long:=0.0; pos.lat:=0.0 END posinval;

PROCEDURE posvalid(pos-:POSITION):BOOLEAN;
BEGIN RETURN (pos.lat<>0.0) OR (pos.long<>0.0) END posvalid;

PROCEDURE radtodeg(p:POSITION):POSITION;
BEGIN
  p.lat:=p.lat*(1.0/RAD);
  p.long:=p.long*(1.0/RAD);
  RETURN p
END radtodeg;


PROCEDURE radmod(w:REAL):REAL;
BEGIN
  WHILE w>PI DO w:=w-2*PI END;
  WHILE w<-PI DO w:=w+2*PI END;
  RETURN w
END radmod;


PROCEDURE wgs84s(lat, long, nn:REAL; VAR x,y,z:REAL);
VAR h, c:REAL;
BEGIN
  h:=nn+EARTH*1000.0;
  z:=h*sin(lat);
  c:=cos(lat);
  y:=h*sin(long)*c;
  x:=h*cos(long)*c;
END wgs84s;

PROCEDURE wgs84sl(lat, long, nn:LONGREAL; VAR x,y,z:LONGREAL);
VAR h, c:LONGREAL;
BEGIN
  h:=nn+EARTH*1000.0;
  z:=h*sin(lat);
  c:=cos(lat);
  y:=h*sin(long)*c;
  x:=h*cos(long)*c;
END wgs84sl;

PROCEDURE wgs84r(x,y,z:REAL; VAR lat, long, heig:REAL);
VAR h:REAL;
BEGIN
  h:=x*x + y*y;
  IF ABS(x)>ABS(y) THEN
    long:=osi.arctan(y/x);
    IF x<0.0 THEN
      IF y>0.0 THEN long:=PI+long ELSE long:=long-PI END;
    END;
  ELSE
    long:=PI*0.5-osi.arctan(x/y);
    IF y<0.0 THEN long:=long-PI END;
  END;
  lat:=osi.arctan(z/sqrt(h));
  heig:=sqrt(h + z*z) - EARTH*1000.0;
END wgs84r;

PROCEDURE wgs84rl(x,y,z:LONGREAL; VAR lat, long:LONGREAL; VAR heig:LONGREAL);
VAR h:LONGREAL;
BEGIN
  h:=x*x + y*y;
  IF ABS(x)>ABS(y) THEN
    long:=atan(y/x);
    IF x<0.0 THEN
      IF y>0.0 THEN long:=PI+long ELSE long:=long-PI END;
    END;
  ELSE
    long:=PI*0.5-atan(x/y);
    IF y<0.0 THEN long:=long-PI END;
  END;
  lat:=atan(z/sqrt(h));
  heig:=sqrt(h + z*z) - EARTH*1000.0;
END wgs84rl;

PROCEDURE distnear(a-,b-:POSITIONL):REAL;
VAR x,y:REAL;
BEGIN
  x:=(a.long-b.long)*coslat;
  y:=a.lat -b.lat;
  RETURN (EARTH*1000.0)*osi.sqrt(x*x+y*y)
END distnear;

PROCEDURE readsymbol(VAR img:SYMBOL; VAR x,y:CARDINAL; fn-:ARRAY OF CHAR);
CONST
      BYTESPERPIX=4;

TYPE 
     ROWS=ARRAY[0..MAXSYMSIZE-1] OF PIX8A;
     pROWS=POINTER TO ROWS;

VAR rows:ARRAY[0..MAXSYMSIZE-1] OF pROWS;
    maxx, maxy, maxxbyte, res, i:INTEGER;
BEGIN
  maxx:=MAXSYMSIZE;
  maxy:=MAXSYMSIZE;
  maxxbyte:=maxx*BYTESPERPIX;                (* maxxbyte = maxx*4: switch on alpha channel *)
  FOR i:=0 TO MAXSYMSIZE-1 DO rows[i]:=ADR(img[i][0]) END; 
  res:=pngread.readpng(fn, rows, maxx, maxy, maxxbyte);
  IF res=0 THEN
    x:=maxx;
    y:=maxy;
  ELSE
    Error("icon image not readable");
--    x:=0;
--    y:=0;
  END;
--WrInt(res, 5); WrInt(maxx, 5); WrInt(maxy, 5); WrStrLn("");
END readsymbol;



PROCEDURE makegammatab(min, max:CARDINAL);
VAR c, span:CARDINAL;
    g:REAL;
BEGIN
  IF max<=min THEN min:=0; max:=HIGH(gammatab) END;
  span:=max-min;
  g:=1.0/igamma;
  FOR c:=0 TO HIGH(gammatab) DO
    IF c>min THEN
      IF c<=max THEN gammatab[c]:=CHR(TRUNC(exp(log(FLOAT(c-min)/FLOAT(span))*g)*255.5));
      ELSE gammatab[c]:=CHR(255) END;
    ELSE gammatab[c]:=CHR(0) END;
  END;
END makegammatab;

PROCEDURE histogram(img:pIMAGE; nosat:BOOLEAN; VAR min, max:CARDINAL);
VAR x, y, pc, pixc, iw, ib, mid, maxe:CARDINAL;
    fmin, mul, fr, fg, fb, w, rw, ilog, llog:REAL;
    h,hh:ARRAY[0..MAXLUM] OF CARDINAL;
--    hr:ARRAY[0..MAXLUM] OF REAL;
BEGIN
  FILL(ADR(h), 0C, SIZE(h));
  FILL(ADR(hh), 0C, SIZE(hh));
  FOR y:=0 TO HIGH(img^[0]) DO
    FOR x:=0 TO HIGH(img^) DO
      WITH img^[x][y] DO
--        pc:=(VAL(CARDINAL, r)*3 + VAL(CARDINAL, g)*5 + VAL(CARDINAL, b)*2) DIV 10; (* PAL luma *)
        pc:=(VAL(CARDINAL, r)*3 + VAL(CARDINAL, g)*5 + VAL(CARDINAL, b)*1) DIV 9;  (* snow shadow blue, dust blue dimms *)
 
        IF nosat THEN
          IF r>pc THEN pc:=r END;                      (* saturated single colours *)
          IF g>pc THEN pc:=g END;
          IF b>pc THEN pc:=b END;
        END;
        IF pc>MAXLUM THEN pc:=MAXLUM END;
        INC(hh[pc]);                                 (* whole image *)
        IF NOT ODD(b) THEN INC(h[pc]) END;           (* only earth no heaven or mirrored heaven *)
      END;
    END;
  END;
  pixc:=(HIGH(img^)+1)*(HIGH(img^[0])+1);
  ib:=truncc(FLOAT(pixc)*ABS(ignorblack)*0.01);
  min:=1;                                                           (* dont use unfilled black stripes *)
  WHILE (ib>=h[min]) & (min<MAXLUM) DO DEC(ib, h[min]); INC(min) END;
  iw:=truncc(FLOAT(pixc)*ignorwhite*0.01);
  max:=MAXLUM;
  WHILE (iw>=hh[max]) & (max>0) DO DEC(iw, hh[max]); DEC(max) END;
  iw:=TRUNC(FLOAT(pixc)*ignorwhite*0.01);
  maxe:=MAXLUM;
  WHILE (iw>=h[maxe]) & (maxe>0) DO DEC(iw, h[maxe]); DEC(maxe) END;

  maxe:=imax(maxe, whitelim);
  min:=imin(min, blacklim);
  max:=imax(max, whitelim);  (* with heaven *)

--- check contrast for auto filmlinearity

  llog:=logfilm;
  IF logfilm<0.0 THEN   (* contrast auto *)
    llog:=200000.0;
    WHILE (llog>100.0) & (llog*osi.ln(1.0+FLOAT(max)/llog)>FLOAT(maxe)) DO llog:=llog*0.95 END;
  END;

  IF llog>0.0 THEN
    ilog:=1.0/llog;
    max:=TRUNC(llog*osi.ln(1.0+FLOAT(max)*ilog));  
  END;

  IF ignorblack<0.0 THEN min:=truncc(FLOAT(max)*ABS(ignorblack)) END;
  IF min>=max THEN max:=min+1 END; 

  fmin:=FLOAT(min);
  mul:=MAXLUM/(FLOAT(max)-fmin);
IF verb THEN WrStr("min..allmax, earthmax, log: "); WrInt(min, 6);WrInt(max, 6); WrInt(maxe, 6);WrFixed(llog, 1,10); WrStrLn(""); END;

  FOR y:=0 TO HIGH(img^[0]) DO
    FOR x:=0 TO HIGH(img^) DO
      WITH img^[x][y] DO
        fr:=FLOAT(r);
        fg:=FLOAT(g);
        fb:=FLOAT(b);

        w:=fr*0.3 + fg*0.5 + fb*0.2 + 0.1;
        rw:=1.0/w;
        fr:=fr*rw;
        fg:=fg*rw;
        fb:=fb*rw;

        w:=w-fmin;
        IF w<0.1 THEN w:=0.1 END;
        IF llog<>0.0 THEN w:=llog*osi.ln(1.0+w*ilog) END;
        w:=mul*w;

        fr:=fr*w;
        IF fr>MAXLUM-1 THEN fr:=FLOAT(MAXLUM-1) END;
        fg:=fg*w;
        IF fg>MAXLUM-1 THEN fg:=FLOAT(MAXLUM-1) END;
        fb:=fb*w;
        IF fb>MAXLUM-1 THEN fb:=FLOAT(MAXLUM-1) END;
        r:=TRUNC(fr);
        g:=TRUNC(fg);
        b:=TRUNC(fb);


  
(*
        r:=TRUNC(w*fr);
        g:=TRUNC(w*fg);
        b:=TRUNC(w*fb);
        IF r>MAXLUM THEN r:=MAXLUM END;
        IF g>MAXLUM THEN g:=MAXLUM END;
        IF b>MAXLUM THEN b:=MAXLUM END;
*)
      END;
    END;
  END;

  min:=0;
  max:=MAXLUM;

END histogram;


PROCEDURE wrpng():INTEGER;

  PROCEDURE pngc(c:CARD16):CARD8;
  BEGIN IF c<=HIGH(gammatab) THEN RETURN ORD(gammatab[c]) ELSE RETURN 255 END END pngc;

VAR pngimg:pngwrite.PNGPIXMAP;
    x, y, ret:INTEGER;
--    min, max:CARDINAL;
BEGIN
  ret:=-2;
  ALLOCATE(pngimg.image, xsize*ysize*3);
  IF pngimg.image<>NIL THEN
(*
    histogram(image, min, max);
    makegammatab(min, max);
*)
    FOR y:=0 TO ysize-1 DO
      FOR x:=0 TO xsize-1 DO
        WITH image^[x][ysize-1-y] DO
          WITH pngimg.image^[x+y*xsize] DO red:=pngc(r); green:=pngc(g); blue:=pngc(b) END;
        END;
      END;
    END;
    pngimg.width:=xsize;
    pngimg.height:=ysize;
    ret:=pngwrite.writepng(imagefn, pngimg);
    DEALLOCATE(pngimg.image, xsize*ysize*3);
  ELSE Werr("png write out of memory"+LF) END;
  RETURN ret
END wrpng;

(*
------- moon

CONST E=RAD*23.4397; (* obliquity of the Earth *)


PROCEDURE fdeclination(l, b:LONGREAL):LONGREAL;
BEGIN RETURN asin(sin(b)*cos(E) + cos(b)*sin(E)*sin(l)) END fdeclination;

PROCEDURE frightascension(l, b:LONGREAL):LONGREAL;
BEGIN RETURN atang2(sin(l)*cos(E)-tan(b)*sin(E), cos(l)) END frightascension;

PROCEDURE siderealtime(d, lw:LONGREAL):LONGREAL;
BEGIN RETURN RAD*(280.16 + 360.9856235*d) - lw END siderealtime;

PROCEDURE faltitude(h, phi, dec:REAL):REAL;
BEGIN RETURN asin(sin(phi)*sin(dec) + cos(phi)*cos(dec)*cos(h)) END faltitude;

PROCEDURE fazimuth(h, phi, dec:REAL):REAL;
BEGIN RETURN atang2(sin(h), cos(h)*sin(phi) - tan(dec)*cos(phi)) END fazimuth;

PROCEDURE mooncoords(d:LONGREAL; VAR declination, distance, rightascension:REAL);
VAR el, m, f, l, b:LONGREAL;
BEGIN
  el:=(218.316 + 13.176396*d)*RAD;
  m:= (134.963 + 13.064993*d)*RAD;
  f:= (93.272  + 13.229350*d)*RAD;

  l:=el + 6.289*RAD*sin(m);         (* latitude *)
  b:=5.128*RAD*sin(f);              (* longitude *)
  distance:=385001 - 20905*cos(m);        (* distance to the moon in km *)
  declination:=fdeclination(l, b);
  rightascension:=frightascension(l, b);
END mooncoords;

PROCEDURE moonposition(lat, long:LONGREAL; time:CARDINAL; VAR az, el, distance:REAL);
VAR lw, phi, declination, bigh, rightascension, h:REAL;
    d:LONGREAL;
s:ARRAY[0..30] OF CHAR;
BEGIN
  lw:=-long;
  phi:=lat;
  d:=(LFLOAT(time) - LFLOAT(86400*((2000-1970)*1461 DIV 4)))/86400.0 - 0.5 (*49925*);
WrFixed(d, 5, 15);WrStrLn("=jd");
  mooncoords(d, declination, distance, rightascension);
  bigh:=siderealtime(d, lw) - rightascension;
  h:=faltitude(bigh, phi, declination);

  (* Altitude correction for refraction *)
  h:=h + (RAD*0.017)/tan(h + (RAD*10.26)/(h + RAD*5.10));

  az:=fazimuth(bigh, phi, declination);
  el:=h;
END moonposition;

PROCEDURE time2moon(t:CARDINAL; mypos:POSITIONL; VAR moonaz, moonel:REAL);
BEGIN
  moonposition(mypos.lat, mypos.long, t, moonaz, moonel, moondist);
  moonaz:=moonaz+PI;
  IF moonaz>2*PI THEN moonaz:=moonaz-2*PI END;
WrStr("moon az el:"); WrFixed(moonaz/RAD, 2, 7); WrFixed(moonel/RAD, 2, 7); WrFixed(moondist, 0, 12);WrStrLn("");
END time2moon;

------- moon
*)

------ sun

PROCEDURE wrsun(sunaz, sunel:REAL);
BEGIN
  IF verb THEN
    WrFixed(sunaz, 2, 7); WrStrLn(" sunaz");
    WrFixed(sunel, 2, 7); WrStrLn(" sunel");
  END;
END wrsun;


PROCEDURE time2sun(date, time:ARRAY OF CHAR; mypos:POSITIONL; VAR sunaz, sunel:REAL):BOOLEAN;
VAR t, day:TIME;

  PROCEDURE dig(c:CHAR; VAR err:BOOLEAN):CARDINAL;
  BEGIN
    IF (c<"0") OR (c>"9") THEN err:=TRUE; RETURN 0 ELSE RETURN ORD(c)-ORD("0") END;
  END dig;

  PROCEDURE StrToDate(s:ARRAY OF CHAR; VAR time:TIME):BOOLEAN; (* <yyyy mm dd> *)
  VAR y, m, d, k:CARDINAL;
      err:BOOLEAN;

  BEGIN
    err:=FALSE;
    y:=dig(s[0], err)*1000 + dig(s[1], err)*100 + dig(s[2], err)*10 + dig(s[3], err);
    IF (y<1970) OR (y>2100) THEN RETURN FALSE END;

    time:=(y-1970)*365 + (y-1969) DIV 4;

    m:=dig(s[4], err)*10 + dig(s[5], err);
    d:=dig(s[6], err)*10 + dig(s[7], err);
    IF (d<1) OR (d>31) THEN RETURN FALSE END;

    k:=ORD(y MOD 4=0);
    INC(time, d-1);
    IF m>2 THEN INC(time, k) END;
    CASE m OF
      1:
     |2: IF d>28+k THEN err:=TRUE END;
         INC(time, 31);
     |3: INC(time, 31+28);
     |4: IF d>30 THEN err:=TRUE END;
         INC(time, 31+28+31);
     |5: INC(time, 31+28+31+30);
     |6: IF d>30 THEN err:=TRUE END;
         INC(time, 31+28+31+30+31);
     |7: INC(time, 31+28+31+30+31+30);
     |8: INC(time, 31+28+31+30+31+30+31);
     |9: IF d>30 THEN err:=TRUE END;
         INC(time, 31+28+31+30+31+30+31+31);
     |10:INC(time, 31+28+31+30+31+30+31+31+30);
     |11:IF d>30 THEN err:=TRUE END;
         INC(time, 31+28+31+30+31+30+31+31+30+31);
     |12:INC(time, 31+28+31+30+31+30+31+31+30+31+30);
     ELSE err:=TRUE;
    END;
    RETURN NOT err
  END StrToDate;

  PROCEDURE StrToDaytime(s:ARRAY OF CHAR; VAR time:TIME):BOOLEAN; (* <yyyy mm dd> *)
  VAR h, mi, se:CARDINAL;
      err:BOOLEAN;

  BEGIN
    err:=FALSE;
    h:=dig(s[0], err)*10 + dig(s[1], err);
    IF h>24 THEN RETURN FALSE END;

    mi:=dig(s[2], err)*10 + dig(s[3], err);
    IF mi>=60 THEN RETURN FALSE END;

    se:=dig(s[4], err)*10 + dig(s[5], err);
    IF se>=60 THEN RETURN FALSE END;

    time:=h*3600 + mi*60 + se;
    RETURN NOT err
  END StrToDaytime;

  PROCEDURE arcsin(x:REAL):REAL;
  BEGIN RETURN osi.arctan(x/sqrt(1.0-x*x)) END arcsin;


  PROCEDURE modrad(w:REAL):REAL;
  BEGIN RETURN (w-FLOAT(TRUNC(w/360.0)*360)) END modrad;

  PROCEDURE sunpos(tt:TIME; lat, long:REAL; VAR azi, ele:REAL);
  VAR l,g,r,aa,e,alpha,c,ga,t0,n,fhg,fg,fl,ta:REAL;
  BEGIN
    n:=FLOAT(tt - 86400*((2000-1970)*1461 DIV 4))/86400.0 - 0.5;
    l:=280.46 + 0.9856474*n;
    l:=modrad(l);

    g:=357.528 + 0.9856003*n;
    g:=modrad(g)*RAD;

    aa:=l + 1.915*sin(g) + 0.02*sin(g*2.0);

    e:=23.439 - 0.0000004*n;

    c:=cos(aa*RAD);
    alpha:=osi.arctan(cos(e*RAD)*sin(aa*RAD)/c)/RAD;
    IF c<0.0 THEN alpha:=alpha + 180.0 END;
    ga:=arcsin(sin(e*RAD)*sin(aa*RAD))/RAD;

    t0:=(VAL(REAL,VAL(INTEGER,n))+0.5)/36525;

    fhg:=6.697376 + 2400.05134*t0 + 1.002738*FLOAT(tt MOD 86400)/3600.0;
    fhg:=fhg - FLOAT(TRUNC(fhg/24.0))*24.0;
    fg:=fhg*15.0;
    fl:=fg+long;
    ta:=fl-alpha;
    c:=cos(ta*RAD)*sin(lat*RAD) - tan(ga*RAD)*cos(lat*RAD);
    azi:=osi.arctan(sin(ta*RAD)/c)/RAD + 360.0 + 180.0;
    IF c<0.0 THEN azi:=azi + 180.0 END;
    azi:=azi-FLOAT(TRUNC(azi/360.0))*360.0;

    ele:=arcsin(cos(ga*RAD)*cos(ta*RAD)*cos(lat*RAD) + sin(ga*RAD)*sin(lat*RAD))/RAD;
    r:=0.017/tan((ele+10.3/(ele+5.11))*RAD);
    ele:=ele+r;
  END sunpos;

  PROCEDURE Refraction(el:REAL):REAL;   (* in rad *)
  BEGIN
    RETURN (PI*0.5-el-arcsin(6366.0/(6366.0+12.0)*cos(el)))/RAD*0.15
  END Refraction;

BEGIN 
  IF StrToDate(date, day) & StrToDaytime(time, t) THEN
    sunpos(t+day*86400, mypos.lat/RAD, mypos.long/RAD, sunaz, sunel);

(*
    sunaz:=sunaz*RAD;
    sunel:=sunel*RAD;

    day:=25568+t DIV (60*60*24);
    day:=1 + day*4 MOD 1461 DIV 4;
 
    AzimutHoehe(mypos.lat, day, FLOAT(t MOD 86400)/3600.0 + mypos.long/(RAD*15.0) - 12.0+zeitgleichung(day), sunaz, sunel);
*)
--WrFixed(sunel,5,10);WrFixed(Refraction(sunel*RAD),5,10);WrStrLn("=sunrefrac");
    sunel:=sunel + Refraction(sunel*RAD);
    wrsun(sunaz, sunel);

--    time2moon(t+day*86400, mypos, moonazim, moonelev);

    RETURN TRUE
 
  END;
  RETURN FALSE;
END time2sun; 

PROCEDURE autosun(VAR az, el:REAL; pos:POSITIONL; pan:REAL);
BEGIN
  IF (pos.lat>0.0)=(pan<180.0) THEN az:=90.0 ELSE az:=-90.0 END;
  az:=pan+az;
  IF az<0.0 THEN az:=az+360.0 ELSIF az>360.0 THEN az:=az-360.0 END;
  el:=15+0.58*(90.0-ABS(pos.lat)/RAD)*ABS(cos(az*RAD));
--WrFixed(ABS(pos.lat)/RAD, 3, 9); WrStrLn("=r");
END autosun;

------- sun

PROCEDURE Parms;
VAR 
    altOG,
    err:BOOLEAN;
    h, iconfn, iconwoodfn, iconhint:ARRAY[0..1023] OF CHAR;
    pan, dist, rn, sunazim, sunelev:REAL;
    label, pri:CARDINAL;
    p:INTEGER;
    ppoi:pPOI;
    ptrk:pTRACK;
    col:INTEGER;

PROCEDURE getpos(VAR p:POSITIONL; e-:ARRAY OF CHAR);
VAR ps:POSITION;
BEGIN
  NextArg(h);
  loctopos(ps, h);
  IF posvalid(ps) THEN p.lat:=ps.lat; p.long:=ps.long
  ELSE
    IF NOT StrToFixL(p.lat, h) OR (ABS(p.lat)>=90) THEN Error(e) END;
    NextArg(h);
    IF NOT StrToFixL(p.long, h) OR (ABS(p.long)>180) THEN Error(e) END;
    p.lat :=p.lat*RAD;
    p.long:=p.long*RAD;
  END; 
END getpos;

BEGIN
  err:=FALSE;
  label:=0;
  iconfn:="";
  iconwoodfn:="";
  iconhint:="";
  pan:=0.0;
  dist:=0.0;
  dateh:="";
  timeh:="";
  altOG:=FALSE;
  millisec:=0;
  pri:=0;
  LOOP
    NextArg(h);
    IF h[0]=0C THEN EXIT END;
    IF (h[0]="-") & (h[1]<>0C) & (h[2]=0C) THEN

     IF h[1]="i" THEN
        NextArg(imagefn);
        IF (imagefn[0]=0C) OR (imagefn[0]="-") THEN Error("-i <imagefilename>") END;
 
     ELSIF h[1]="P" THEN
        ALLOCATE(ppoi, SIZE(ppoi^));
        IF ppoi=NIL THEN Error("poifile, out of memory") END;

        ppoi^.r:=255;
        ppoi^.g:=255;
        ppoi^.b:=255;
        ppoi^.alpha:=255;
        NextArg(h);
        ppoi^.prio:=pri;                         (* prefer first pois if no place for text *)
        INC(pri);
        IF StrToInt(h, col) THEN
          ppoi^.r:=col;
          NextArg(h);
          IF StrToInt(h, col) THEN
            ppoi^.g:=col;
            NextArg(h);
            IF StrToInt(h, col) THEN
              ppoi^.b:=col;
              NextArg(h);
              IF StrToInt(h, col) THEN
                ppoi^.alpha:=col;
                NextArg(h);
              END;
            END;
          END;
        END; 

        Assign(ppoi^.poifn, h);
        IF (ppoi^.poifn[0]=0C) OR (ppoi^.poifn[0]="-") THEN Error("-P <POI filename>") END;
        Assign(ppoi^.iconfn, iconfn);
        Assign(ppoi^.iconwoodfn, iconwoodfn);
        Assign(ppoi^.iconhint, iconhint);
        ppoi^.altIsOG:=altOG;
        altOG:=FALSE;
        ppoi^.next:=poifiles;
        poifiles:=ppoi;        

     ELSIF h[1]="E" THEN
        ALLOCATE(ptrk, SIZE(ptrk^));
        IF ptrk=NIL THEN Error("trackfile, out of memory") END;

        ptrk^.r:=200;
        ptrk^.g:=100;
        ptrk^.b:=50;
        ptrk^.width:=2;
        ptrk^.altmode:=0;
        NextArg(h);
        IF StrToInt(h, col) THEN
          ptrk^.r:=col;
          NextArg(h);
          IF StrToInt(h, col) THEN
            ptrk^.g:=col;
            NextArg(h);
            IF StrToInt(h, col) THEN
              ptrk^.b:=col;
              NextArg(h);
              IF StrToInt(h, col) THEN
                ptrk^.width:=col;
                NextArg(h);
                IF StrToInt(h, col) THEN
                  ptrk^.altmode:=col;
                  NextArg(h);
                END;
              END;
            END;
          END;
        END; 

        Assign(ptrk^.fn,h);    (* for -v *)
        ptrk^.fd:=OpenRead(h);
        IF ptrk^.fd<0 THEN Error("-E trackfile not readable") END;
        IF (h[0]=0C) OR (h[0]="-") THEN Error("-E <trackfilename>") END;
        ptrk^.next:=trackfiles;
        trackfiles:=ptrk;        

     ELSIF h[1]="c" THEN
        NextArg(csvfn);
        IF (csvfn[0]=0C) OR (csvfn[0]="-") THEN Error("-c <csv filename>") END;

     ELSIF h[1]="C" THEN
        NextArg(h);
        IF NOT StrToInt(h, col) OR (col<0) THEN Error("-C <number>") END;
        rastercsv:=col;

     ELSIF h[1]="I" THEN
        NextArg(iconfn);
        IF (iconfn[0]=0C) OR (iconfn[0]="-") THEN Error("-I <icon filename> <icon no wood filename> (png)") END;
        NextArg(iconwoodfn);
        IF iconwoodfn[0]="-" THEN Error("-I <icon filename> <icon no wood filename> (png)") END;

     ELSIF h[1]="J" THEN
        NextArg(iconhint);
        IF (iconhint[0]=0C) OR (iconhint[0]="-") THEN Error("-J <icon hint>") END;
        IF iconhint[0]="," THEN iconhint[0]:=0C END;

     ELSIF h[1]="v" THEN verb:=TRUE;
     ELSIF h[1]="O" THEN altOG:=TRUE;
     ELSIF h[1]="Z" THEN defaultocean:=TRUE;
     ELSIF h[1]="j" THEN dropshadows:=TRUE;
     ELSIF h[1]="u" THEN
       NextArg(h);
       IF NOT StrToFix(urbanoff, h) THEN Error("-u <pixel per meter>") END;

     ELSIF h[1]="f" THEN
       NextArg(h);
       IF h[0]="0" THEN flatscreen:=eFLAT ELSIF h[0]="1" THEN flatscreen:=eCAMERA ELSE Error("-f <n>") END;

     ELSIF h[1]="x" THEN
        NextArg(h);
        IF NOT StrToInt(h, xsize) OR (xsize<20) OR (xsize>=32000) THEN Error("-x <size>") END;

     ELSIF h[1]="y" THEN
        NextArg(h);
        IF NOT StrToInt(h, ysize) OR (ysize<20) OR (ysize>=8192) THEN Error("-y <size>") END;

     ELSIF h[1]="M" THEN
        NextArg(h);
        IF NOT StrToCard(h, srtmmaxmem) THEN Error("-M <bytes>") END;

     ELSIF h[1]="p" THEN
        NextArg(srtmdir);
        IF (srtmdir[0]=0C) OR (srtmdir[0]="-") THEN Error("-p <srtm-folder-path>") END;

     ELSIF h[1]="a" THEN getpos(posa, "-a <lat> <long> or <locator>"); 

     ELSIF h[1]="b" THEN getpos(posb, "-b <lat> <long> or <locator>");
       
     ELSIF h[1]="A" THEN
        NextArg(h);
        IF NOT StrToFix(camog, h) OR (camog<0) OR (camog>500000) THEN Error("-A <meter> (0..500000)") END;

     ELSIF h[1]="d" THEN
        NextArg(h);
        IF NOT StrToFix(pan, h) THEN Error("-d <deg> <km>") END;
        NextArg(h);
        IF NOT StrToFix(dist, h) THEN Error("-d <deg> <km>") END;

     ELSIF h[1]="S" THEN
        NextArg(dateh);
        NextArg(timeh);

     ELSIF h[1]="V" THEN
       NextArg(h);
       IF NOT StrToCard(h, poilabel) THEN Error("-V <n>") END;
       INC(poilabel);

     ELSIF h[1]="W" THEN
        flatwater:=TRUE;
        NextArg(h);
        IF NOT StrToFix(waveamp, h) OR (waveamp<0.0) OR (waveamp>10.0) THEN Error("-W <wavehigth> <wavelength> <flatness>") END;
        waveamp:=waveamp*0.01;
        NextArg(h);
        IF NOT StrToFix(wavelength, h) THEN Error("-W <wavehigth> <wavelength> <flatness>") END;
        NextArg(h);
        IF NOT StrToFix(waterslant, h) THEN Error("-W <wavehigth> <wavelength> <flatness>") END;


     ELSIF h[1]="m" THEN
        NextArg(h);
        IF NOT StrToFix(rn, h) THEN Error("-m <maxmountainhight>"); END;
        maxmountain:=rn;

     ELSIF h[1]="N" THEN
        NextArg(h);
        IF NOT StrToFix(abscolmin, h) THEN
          Assign(colfn, h);
          NextArg(h);
          IF NOT StrToFix(abscolmin, h) THEN Error("-N [<coloufile>] <minalt> <maxalt> <step>"); END;
        END;
        NextArg(h);
        IF NOT StrToFix(abscolmax, h) THEN Error("-N [<coloufile>] <minalt> <maxalt> <step>"); END;
        NextArg(h);
        IF NOT StrToCard(h, altcolstep) THEN Error("-N [<coloufile>] <minalt> <maxalt> <step>"); END;
        alticolour:=TRUE;

     ELSIF h[1]="t" THEN
        NextArg(h);
        IF NOT StrToFix(treesize, h) THEN Error("-t <treesize> <botanicslant> <treedist>"); END;
        NextArg(h);
        IF NOT StrToFix(botanicslant, h) THEN Error("-t <treesize> <botanicslant> <treedist>"); END;
        NextArg(h);
        IF NOT StrToFix(treerastersize, h) OR (treerastersize<1.0) OR (treerastersize>100.0) THEN Error("-t <treesize> <botanicslant> <treedist> [1..100m]") END;

     ELSIF h[1]="T" THEN
        NextArg(h);
        IF NOT StrToFix(zerotreealt, h) THEN Error("-T <no-tree-alt> <full-tree-alt>"); END;
        NextArg(h);
        IF NOT StrToFix(fulltreealt, h) THEN Error("-T <no-tree-alt> <full-tree-alt>"); END;
        IF zerotreealt<=fulltreealt THEN Error("-T no_tree_alt <= full_tree_alt ?"); END;

     ELSIF h[1]="l" THEN
        NextArg(h);
        IF NOT StrToFix(ignorblack, h) THEN Error("-l <black%> <white%>"); END;
        NextArg(h);
        IF NOT StrToFix(ignorwhite, h) THEN Error("-l <black%> <white%>"); END;

     ELSIF h[1]="D" THEN
        NextArg(h);
        IF NOT StrToFix(desertcol.r, h) THEN Error("-D <r> <g> <b> <granular>"); END;
        NextArg(h);
        IF NOT StrToFix(desertcol.g, h) THEN Error("-D <r> <g> <b> <granular>"); END;
        NextArg(h);
        IF NOT StrToFix(desertcol.b, h) THEN Error("-D <r> <g> <b> <granular>"); END;
        NextArg(h);
        IF NOT StrToFix(sandgranularity, h) THEN Error("-D <r> <g> <b> <granular>"); END;
(*
        desertcol.r:=desertcol.r*16.0;
        desertcol.g:=desertcol.g*16.0;
        desertcol.b:=desertcol.b*16.0;
*)
     ELSIF h[1]="o" THEN
        NextArg(h);
        IF NOT StrToFix(minpoialt, h) OR (minpoialt<-20000.0) OR (minpoialt>20000.0) THEN Error("-o <meter> (0..20000)") END;

     ELSIF h[1]="K" THEN
        NextArg(h);
        IF NOT StrToFix(logfilm, h) THEN Error("-K <filmlinearity>") END;

     ELSIF h[1]="G" THEN
        NextArg(h);
        IF NOT StrToFix(snowalt, h) THEN Error("-G <meter> <meter>") END;
        NextArg(h);
        IF NOT StrToFix(slowaltslant, h) THEN Error("-G <meter> <meter>") END;
        IF slowaltslant<1.0 THEN slowaltslant:=1.0 END;

     ELSIF h[1]="r" THEN
        NextArg(h);
        IF NOT StrToFix(refraction, h) OR (refraction<0.0) OR (refraction>1.0)
        THEN Error("-r <refraction> [0.0..1.0]") END;

     ELSIF h[1]="e"  THEN
        NextArg(h);
        IF NOT StrToFix(elev, h) OR (elev<-90.0) OR (elev>90.0)
        THEN Error("-e <elevation> [-90.0..90.0]") END;

     ELSIF h[1]="w"  THEN
        NextArg(h);
        IF NOT StrToFix(angle, h) OR (angle<0.1) OR (angle>360.0)
        THEN Error("-w <degrees> [0.1..360.0]") END;

     ELSIF h[1]="L"  THEN
        NextArg(h);
        IF NOT StrToCard(h, blacklim) THEN Error("-L <blacklim> <whitelim>"); END;
        NextArg(h);
        IF NOT StrToCard(h, whitelim) THEN Error("-L <blacklim> <whitelim>"); END;

     ELSIF h[1]="Q"  THEN
        NextArg(h);
        IF NOT StrToFix(videolat, h) OR (videolat<-89.0) OR (videolat>89.0)
        THEN Error("-Q <degrees> [-89..89]") END;

     ELSIF h[1]="z"  THEN
        NextArg(h);
        IF NOT StrToFix(yzoom, h) OR (yzoom<0.001) OR (yzoom>10.0)
        THEN Error("-z <factor> [0.0..10.0]") END;

     ELSIF h[1]="s" THEN
        NextArg(h);
        IF NOT StrToCard(h, poisize) OR (poisize>20) THEN Error("-s <size> 1..20") END;

     ELSIF h[1]="g" THEN
        NextArg(h);
        IF NOT StrToFix(igamma, h) OR (igamma<0.1) OR (igamma>10.0) THEN Error("-g <gamma> [0.1..10] (1.0)") END;

     ELSIF h[1]="H" THEN
        NextArg(h);
        IF NOT StrToFix(heavencol.r,h) OR (heavencol.r>=256.0) THEN Error("-H <r> <g> <b> 0..255") END;
        NextArg(h);
        IF NOT StrToFix(heavencol.g,h) OR (heavencol.g>=256.0) THEN Error("-H <r> <g> <b> 0..255") END;
        NextArg(h);
        IF NOT StrToFix(heavencol.b,h) OR (heavencol.b>=256.0) THEN Error("-H <r> <g> <b> 0..255") END;
        heavencol.r:=heavencol.r*16.0;
        heavencol.g:=heavencol.g*16.0;
        heavencol.b:=heavencol.b*16.0;

     ELSIF h[1]="R" THEN
        NextArg(h);
        IF NOT StrToFix(watercol.r,h) OR (watercol.r>=1024.0) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>") END;
        NextArg(h);
        IF NOT StrToFix(watercol.g,h) OR (watercol.g>=1024.0) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>") END;
        NextArg(h);
        IF NOT StrToFix(watercol.b,h) OR (watercol.b>=1024.0) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>") END;
        NextArg(h);
        IF NOT StrToFix(waterglance,h) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>>") END;
        NextArg(h);
        IF NOT StrToFix(watertransparency,h) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>") END;
        NextArg(h);
        IF NOT StrToFix(surfacedust,h) THEN Error("-R <r> <g> <b> <glance> <transp> <surface>") END;

     ELSIF h[1]="F" THEN
        NextArg(h);
        IF h[0]="1" THEN fonttyp:=6 ELSIF h[0]="2" THEN fonttyp:=8 ELSIF h[0]="3" THEN fonttyp:=10
        ELSE Error("-F <font>[,<font>] (1:6x10 2:8x14 3:10x20)") END;
        IF h[1]="," THEN
        IF h[2]="1" THEN fonttypfoot:=6 ELSIF h[2]="2" THEN fonttypfoot:=8
        ELSIF h[2]="3" THEN fonttypfoot:=10
        ELSE Error("-F <font>[,<font>] (1:6x10 2:8x14 3:10x20)") END;

        END;

     ELSIF h[1]="h" THEN
        WrStrLn("");
        WrStrLn("Panorama");
        WrStrLn(" -A <m>                            Camera over ground or over NN [m] (10)");
        WrStrLn(" -a <lat> <long> | [locator]       Camera position lat long (degrees) or qth locator");
        WrStrLn(" -b <lat> <long> | [locator]       Pan to point lat long (degrees) or qth locator");
        WrStrLn(" -C <N>                            write CSV line for every N'th column");
        WrStrLn(" -c <filename>                     CSV file name");
        WrStrLn(" -D <r> <g> <b> <granularity>      Sand/Rock colour and granularity (with -t 0 0 -u 0)");
        WrStrLn(" -d <deg> <km>                     Relative to position a (alternativ to -b)");
        WrStrLn("                                     km distance sight limit");
        WrStrLn(" -E [<r> <g> <b> <width> <mode>] <trackfile> draw track (lat,long[,altNN]) mode=0: always on ground,");
        WrStrLn("                                     mode=1: use alt if above ground, mode=2: no track if alt under gnd");
        WrStrLn("                                     mode=3: use altOG with minimum -o <m>");
        WrStrLn("                                     blank csv line starts a new track (same colour/mode)");
        WrStrLn(" -e <degrees>                      Camera Elevation (degrees) (0.0)");
        WrStrLn("                                     speed optimizer, surface contrast and dust calculation");
        WrStrLn("                                     are able to handle only small elevations!");
        WrStrLn(" -F <font>[,<font>]                Font Size POI,footline (1) 1: 6x10, 2: 8x14, 3: 10x20");
        WrStrLn(" -f <n>                            0 Flat screen perpendicular, 1 screen normal to sight (no latitude scale)");
        WrStrLn("                                     default courved projection area (needed at >=180deg sight)");
        WrStrLn(" -G <meter> <meter>                Glazier/Snow altitude and fade out (5000 500)");
        WrStrLn(" -g <gamma>                        Image Gamma 0.1..10 (1.0)");
        WrStrLn(" -H <r> <g> <b>                    Heaven colour (50 70 300)");
        WrStrLn(" -h                                this");
        WrStrLn(' -I <filename> <filename> | ""     Symbol Image File Name with, w/o wood, set before -P');
        WrStrLn(" -i <filename>                     Image File Name");
        WrStrLn(" -j                                dropshadows on (test implementation, very slow)");
        WrStrLn(" -J <poi-hint>                     text appended to csv lines, set before -P, off with -J ,");
        WrStrLn(" -K <filmcharacteristic>           compress sunlight/dust to fit in image brightness (4000) 0=linear else auto");
        WrStrLn(" -l <%> <%>                        fit image contrast to saturate % of all black/white pixels (-0.1 1.0)");
        WrStrLn("                                   negative black level use -level of while");
        WrStrLn(" -L <n> <n>                        limit contrast expansion black white (400 1000)");
        WrStrLn(" -M <bytes>                        SRTM-Cache Limit (100000000)");
        WrStrLn(" -m <meter>                        limit searching montain size in viewed area for speed on long distance (8900)");
        WrStrLn(" -N [<colourfile>] <meter> <meter> <meter>  Colourize altitude <from><to><step> m, 0 is autoscale");
        WrStrLn("                                      file: [<a>] <r> <g> <b> ... min to max (0..255)");
        WrStrLn("                                      <a>: w water, t trees, u urban, x ?, y ?");
        WrStrLn(" -O                                POI File altitude is over ground not NN, set before -P");
        WrStrLn(" -o <m>                            min. POI altitude over ground (1)");
        WrStrLn("                                     if given in POI file and higher use this");
        WrStrLn(" -P [<r> <g> <b> [<transp>]] <filename>  optional icon colours (0..255) POI File Name");
        WrStrLn("                                     repeat for more files add -V for labeling");
        WrStrLn(" -p <pathname>                     folder with /srtm1 /srtm3 /srtm30");
        WrStrLn(" -Q <lat>                          for video enter raw latitude for squarish tree raster (default camera lat)");
        WrStrLn("                                     camera altitude now over NN");
        WrStrLn(" -R <r> <g> <b> <glance> <transparency> Water Colour mirroring sight  (350 300 200 1 1)");
        WrStrLn(" -r <refraction>                   0.0(vacuum)..1.0(earth is a disk) (0.13)");
        WrStrLn(" -S <deg> <deg>[+<ms>] | <date> <time>[+<ms>]  Sun azimuth and elevation or utc: 20191231 115959");
        WrStrLn("                                   + append milliseconds relative time for waves in video 20191231 115959+8000");
        WrStrLn(" -s <size>                         internal generated POI symbol size (5)");
        WrStrLn(" -T <no-tree-alt> <full-tree-alt>  2500 1200 m altitude shrink trees from full size to <2m");
        WrStrLn(" -t <meter> <steep> <dist>         Enable vegetation, Tree size (m) lowered linear as given in -T, 0=wood off");
        WrStrLn("                                     steep=0 all rock, steep=5 rock if steeper");
        WrStrLn("                                     median tree distance in wood (8)");
        WrStrLn(" -u <meter>                        Urban area fade-in at pixel per meter, 0 urban off");
        WrStrLn(" -V <n>                            Write POI-Name over Icon, 0 all, 1 as place for text");
        WrStrLn("                                     sortet prior on lefter -P and than nearer objects");
        WrStrLn(" -v                                Say something");
        WrStrLn(" -W <wavehigth> <wavelength> <flatness>  draw Water with waves (0.3) (1.0) (0.0), flatness 0 no filter, 0.99 very flat");
        WrStrLn("                                     -0.99 make water on flat if no tagged srtm");
        WrStrLn(" -w <degrees>                      Camera horizontal sight angle (degrees) (45.0)");
        WrStrLn(" -x <size>                         Image size x (600)");
        WrStrLn(" -y <size>                         Image size y (400)");
        WrStrLn(" -Z                                Ocean if no (valid) srtm data or srtm30 with 0m");
        WrStrLn(" -z <factor>                       Vertical Zoom (for 360deg panorama) (1.0)");
        WrStrLn("");
        HALT

      ELSE err:=TRUE END;
    ELSE err:=TRUE END;
    IF err THEN EXIT END;
  END;
  IF err THEN Werr(">"); Werr(h); Werr("< use -h"+LF); HALT END;

--  IF (csvfn[0]<>0C) & (flatscreen=eCAMERA) THEN Werr("no map csv on not vertical canvas"+NL) END;

  IF altOG & (poifiles<>NIL) & (poifiles^.next=NIL) THEN poifiles^.altIsOG:=altOG END; (* exact 1 poi file so altOG switch may be anywhere *)

  p:=InStr(timeh, "+");
  IF p>0 THEN
    timeh[p]:=0C;
    INC(p);
    WHILE (p<=HIGH(timeh)) & (timeh[p]<>0C) DO
      millisec:=millisec*10 + (ORD(timeh[p])-ORD("0"));
      INC(p);
    END;
  END;
--WrInt(millisec, 10); WrStrLn("");

  IF time2sun(dateh, timeh, posa, sunazim, sunelev) THEN
--    dateh:=""; timeh:="";
    sunaz:=sunazim;
    sunel:=sunelev;

  ELSIF StrToFix(sunazim, dateh) & StrToFix(sunelev, timeh) THEN
    IF ABS(sunelev)>90.0 THEN Error("-S <deg> <deg> or yyyymmdd hhmmss") END;

    dateh:=""; timeh:="";
    sunaz:=sunazim;
    sunel:=sunelev;
  END;

  IF NOT posvalidl(posb) & posvalidl(posa) & (dist>0.01) THEN
    pan:=pan*RAD;
    posb.lat:=posa.lat + cos(pan)*dist/EARTH;
    posb.long:=posa.long + sin(pan)/cos(posa.lat)*dist/EARTH;
--WrFixed(posb.lat/RAD, 5,10); WrFixed(posb.long/RAD, 5,10);
--WrFixed(azimuth(posa, posb), 2,10);WrFixed(distance(posa, posb), 3,10); WrStrLn("");
    IF NOT posvalidl(posb) THEN Error("camera sight point out of database") END;
  END;

  IF sunel<-90.0 THEN
    autosun(sunaz, sunel, posa, azimuthl(posa, posb));
--    IF verb THEN Werr('no "-S",  default sun position generatet'+LF) END;
    wrsun(sunaz, sunel); 
  END;

  treeshrink:=1.0/(zerotreealt-fulltreealt); 

  IF alticolour THEN
    treesize:=0.0;
    IF logfilm<0.0 THEN logfilm:=0.0 END;
  END;
END Parms;

PROCEDURE readcolourfile(VAR col:ALTCOLTAB; fn-:ARRAY OF CHAR);
VAR fd:File;
    fb:ARRAY[0..30000] OF CHAR;
    b:ARRAY[0..100] OF CHAR;
    l, i, p, j, x, tag:INTEGER;
    hc:CCOL;
BEGIN
  FILL(ADR(col), 0C, SIZE(col));
  IF fn[0]=0C THEN RETURN END;

  fd:=OpenRead(fn);
  IF fd<0 THEN Error("Colourtable file not found") END;

  l:=RdBin(fd, fb, SIZE(fb));
  IF l<=0 THEN Error("Colourtable file not readable") END;
  IF l=SIZE(fb) THEN Error("Colourtable too long") END;

  p:=0;
  j:=0;
  tag:=-1;
  LOOP
    i:=0;
    WHILE (p<l) & (fb[p]<=" ") DO INC(p) END;
    WHILE (p<l) & (i<VAL(INTEGER, HIGH(b))) & (fb[p]>" ") DO
      b[i]:=fb[p];
      INC(i);
      INC(p);
    END;
    b[i]:=0C;
    IF col.len>HIGH(col.c) THEN Error("too much Colours") END;

    IF i>0 THEN
      IF (j=0) & (b[0]="w") THEN tag:=0
      ELSIF (j=0) & (b[0]="t") THEN tag:=1
      ELSIF (j=0) & (b[0]="u") THEN tag:=2
      ELSIF (j=0) & (b[0]="s") THEN tag:=3
      ELSIF (j=0) & (b[0]="x") THEN tag:=4
      ELSIF (j=0) & (b[0]="y") THEN tag:=5
      ELSE
        IF NOT StrToInt(b, x) THEN Error("Junc in Colourfile") END;

        IF x<0 THEN x:=0 ELSIF x>255 THEN x:=255 END;
        IF j=0    THEN hc.r:=x; INC(j);
        ELSIF j=1 THEN hc.g:=x; INC(j);
        ELSE 
          hc.b:=x; 
          j:=0;
          IF tag>=0 THEN col.tag[tag]:=hc; tag:=-1 ELSE col.c[col.len]:=hc; INC(col.len) END;
        END;
      END;
    ELSE EXIT END;

  END;
END readcolourfile;


PROCEDURE wrcsv(px, py:INTEGER; posdeg-:POSITION; alti:REAL; text-, iconhint-, ext-:ARRAY OF CHAR);
VAR s,h:ARRAY[0..1000] OF CHAR;
BEGIN
  IF csvfn[0]<>0C THEN
    IF csvfd<0 THEN csvfd:=OpenWrite(csvfn) END;
    IF csvfd<0 THEN Error("cannot write csv-file") END;

    IntToStr(px, 0, s); Append(s,",");
    IntToStr(py, 0, h); Append(s,h);  Append(s,",");
    FixToStr(posdeg.lat , 6,  h); Append(s,h); Append(s,",");
    FixToStr(posdeg.long, 6,  h); Append(s,h); Append(s,",");
    IF (alti>-1000.0) & (alti<10000.0) THEN
      IntToStr(VAL(INTEGER,alti+0.5), 0,  h); Append(s,h);
    END;
    Append(s,',"');
    Append(s, text);
    Append(s,'",');
    Append(s, iconhint);
    Append(s,',"');
    Append(s, ext);  
    Append(s, '"'+NL);
    WrBin(csvfd, s, Length(s));
  END;


-- WITH image^[px][py] DO INC(r, 200); INC(g, 200); END;


END wrcsv;



PROCEDURE xycsv(px, py:CARDINAL; posl:POSITIONL; flush, cont:BOOLEAN);
CONST DDEG=50.0/20000000.0*PI;
      DALT=10.0;
VAR alt, resol, dx, x0, y0, a0, kx, ky, ka, nf:REAL;
    n, i:CARDINAL;
    eov:BOOLEAN;
    pos:POSITION;
BEGIN
  pos.long:=posl.long;
  pos.lat:=posl.lat;
  alt:=getsrtm(pos, 0, resol);
  IF defaultocean & (alt>NOALT) THEN alt:=0.0 END; 
  IF cont THEN
    n:=csvvec.act+1;
    IF (n>2) OR (n>0) & flush THEN
      IF flush OR (n>=HIGH(csvvec.v)) THEN eov:=TRUE 
      ELSE
        nf:=FLOAT(n);
        x0:=csvvec.v[0].p.long;
        y0:=csvvec.v[0].p.lat;
        a0:=csvvec.v[0].a; 
        kx:=(pos.long-x0)/nf;
        ky:=(pos.lat -y0)/nf;
        ka:=(alt-a0)/FLOAT(n);
        eov:=FALSE; 
        dx:=DDEG*cos(pos.lat);      
        i:=n;
        REPEAT 
          DEC(i);
          nf:=FLOAT(i);
          IF (i>=n) OR (ABS(x0 + kx*nf - csvvec.v[i].p.long)>dx)
          OR (ABS(y0 + ky*nf - csvvec.v[i].p.lat )>DDEG)
          OR (ABS(a0 + ka*nf - csvvec.v[i].a)>DALT) THEN eov:=TRUE END;
        UNTIL (i=0) OR eov;
      END;
      IF eov THEN
        wrcsv(px, py-n, radtodeg(csvvec.v[0].p), csvvec.v[0].a, "", "", "");
        IF n>1 THEN
          wrcsv(px, py-1, radtodeg(csvvec.v[n-1].p), csvvec.v[n-1].a, "", "", "");
        END;
        n:=0; 
      END;
    END; 
  ELSE n:=0 END; 
  csvvec.act:=n;
  csvvec.v[n].a:=alt;
  csvvec.v[n].p:=pos;
END xycsv;


PROCEDURE AddShadow(dist, alt:LONGREAL);
VAR s:pDROPSHADOW;
BEGIN
  IF (newshadow<>NIL) & (newshadow^.dist>=dist) THEN      (* trace goes back *) 
    IF newshadow^.alt<alt THEN newshadow^.alt:=alt END;
  ELSE
    IF freeshadow<>NIL THEN
      s:=freeshadow;
      freeshadow:=freeshadow^.next;
    ELSE ALLOCATE(s, SIZE(s^)) END;
    IF s<>NIL THEN
      s^.next:=newshadow;
      newshadow:=s;
      s^.dist:=dist;
      s^.alt:=alt;
    END;
  END;
END AddShadow;

PROCEDURE Shadowmask(xi:INTEGER);
VAR
BEGIN

END Shadowmask;


PROCEDURE colblend(v:REAL; c0-,c1-:COL; VAR c:COL);
VAR vv:REAL;
BEGIN
  IF v<0.0 THEN v:=0.0 ELSIF v>1.0 THEN v:=1.0 END;
  vv:=1.0-v;
  c.r:=c0.r*vv + c1.r*v;
  c.g:=c0.g*vv + c1.g*v;
  c.b:=c0.b*vv + c1.b*v;
END colblend;


PROCEDURE azieletorotvec(az,el:REAL; VAR n:ROTVEC);
BEGIN
  n.sx:=osi.sin(az);
  n.cx:=osi.cos(az);
  n.sy:=osi.sin(el);
  n.cy:=osi.cos(el);
END azieletorotvec;

PROCEDURE azieletovec(az,el:REAL; VAR n:NORMVEC);
VAR r:REAL;
BEGIN
  r:=osi.cos(el);
  n.x:=r*osi.sin(az);
  n.y:=r*osi.cos(az);
  n.z:=osi.sin(el);
END azieletovec;




PROCEDURE normvec(pos-:POSITIONL; pmeta:pMETAINFO; VAR slant, h:REAL; VAR nvec:NORMVEC;
                  VAR water, wood, urban:CARD8);
CONST TESTDIST=10.0;
VAR 
    sqxy,
    resol   :REAL;
    attr    :CARD8;
BEGIN
  pmeta^.withslant:=TRUE;
  pmeta^.aliasattr:=TRUE;
  h:=getsrtmlong(pos.lat, pos.long, 1, TRUE, resol, attr, pmeta);
  IF defaultocean & ((h>NOALT) OR (h=0.0) & (resol>=900.0)) THEN
    h:=0.0;
    pmeta^.attrweights[ATTRWATER]:=1.0;
  END;
  water:=TRUNC(pmeta^.attrweights[ATTRWATER]*100.0);
  IF (treesize>0.0) OR alticolour THEN wood:=TRUNC(pmeta^.attrweights[ATTRWOOD]*100.0) ELSE wood:=0 END;
  urban:=TRUNC(pmeta^.attrweights[ATTRURBAN]*100.0);
--  IF NOT alticolour THEN
--dir:=atang2(meta.slantx, meta.slanty);
--slant:=1.0-osi.arctan(osi.sqrt(sqr(meta.slantx)+sqr(meta.slanty)))*(0.5/PI);
  sqxy:=sqr(pmeta^.slantx) + sqr(pmeta^.slanty);
  slant:=osi.cos(osi.arctan(osi.sqrt(sqxy)));                     (* do it faster *)
--  END;
  nvec.z:=1.0/osi.sqrt(sqxy + 1.0);
  nvec.x:=pmeta^.slantx*nvec.z;
  nvec.y:=pmeta^.slanty*nvec.z;
END normvec;

(*
PROCEDURE normvectree(VAR th:TREEHINT);
BEGIN
  WITH th DO
    IF hitradius<1.0 THEN
      nvec.x:=-hitx; nvec.y:=-hity; nvec.z:=1.0-osi.sqrt(hitradius)*treesize/peakh;
    ELSE nvec.x:=0.0; nvec.y:=0.0; nvec.z:=1.0; END;
  END;
END normvectree;
*)

PROCEDURE SunRefraktion(ele:REAL):REAL;
VAR r, t, p:REAL;
BEGIN
(*
 Naeherungsloesung für die Refraktion fuer ein Objekt bei Hoehe hoehe über mathematischem Horizont
 Refraktion beträgt bei Sonnenaufgang 34 Bogenminuten = 0.56667°
 Falls die Höhe der Sonne nicht genauer als auf 0.5° gewuenscht ist, kann diese Funktion ignoriert werden
*)
  p:=1013.25;       (* Luftdruck der Standard-Atmosphäre in hPa (=mbar) *)
  t:=15.0;          (* Temperatur der Standard-Atmosphäre in gradC *)
  r:=0.0;
  IF ele>=15.0*RAD THEN r:=0.00452*RAD*p/tan(ele)/(273.0+t);  (* ueber 15grad - einfachere Formel *)
  ELSIF ele>-1.0*RAD THEN r:=RAD*p*(0.1594+0.0196*ele+0.00002*ele*ele)/((273.0+t)*(1+0.505*ele+0.0845*ele*ele)) END;
  RETURN r          (* Refraktion in Radians *)
END SunRefraktion;












PROCEDURE suncolour(sunele:REAL; VAR cfar:COL; VAR sunlum:REAL);

VAR s, es:REAL;
BEGIN  
  s:=1.0-sunele*(1.0/(90*RAD));                              (* sunlight colour *)  
  IF s<0.001 THEN s:=0.001 END;
  es:=pow(s,12.0);
  IF es>1.0 THEN es:=1.0/es END;
  cfar.r:=cFARCHROMr-es*(0.1*cFARCHROMr);
  cfar.g:=cFARCHROMg-es*(0.6*cFARCHROMg);
  cfar.b:=cFARCHROMb-es*(1.0*cFARCHROMb);


  sunlum:=cfar.r*0.3+cfar.g*0.5+cfar.b*0.2;                  (* luminance part of sun *)
  IF sunele<0.0 THEN
    s:=1.0+sunele*5.0;
    IF s<0.05 THEN s:=0.05 END;
    sunlum:=sunlum*s;
  END;
END suncolour;


PROCEDURE shadowcolour(shadowval:REAL; suncol-:COL; VAR shadwocol:COL);
BEGIN
  IF shadowval<0.0 THEN shadowval:=0.0 ELSIF shadowval>1.0 THEN shadowval:=1.0 END;
--  IF tree THEN shadowval:=shadowval*0.4+0.3 END;
  shadwocol.r:=suncol.r*cMAXCHROMr*shadowval;
  shadwocol.g:=suncol.g*cMAXCHROMg*shadowval;
  shadwocol.b:=suncol.b*cMAXCHROMb*shadowval;
END shadowcolour;






PROCEDURE irand(x, y:LONGREAL):REAL;                    (* bilinear interpolated random map *)
VAR a,b,c,d,fx,fy:REAL;
    ix,iy,ixx,iyy:CARDINAL;
BEGIN
  x:=ABS(x);
  y:=ABS(y);
  ix:=TRUNC(x);
  fx:=x-LFLOAT(ix);
  iy:=TRUNC(y);
  fy:=y-LFLOAT(iy);

  ix:=ix MOD NOISESIZE;
  iy:=iy MOD NOISESIZE;
  ixx:=(ix+1) MOD NOISESIZE;
  iyy:=(iy+1) MOD NOISESIZE;
  
  a:=noisemap[ix ][iy ];
  b:=noisemap[ixx][iy ];
  c:=noisemap[ix ][iyy];
  d:=noisemap[ixx][iyy];

  a:=a*(1.0-fx) + b*fx;
  b:=c*(1.0-fx) + d*fx;
  RETURN a*(1.0-fy) + b*fy
END irand;


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 qrand(x, y:LONGREAL):REAL;                (* biqubic interpolated random map *)
VAR fx,fy:REAL;
    ix,iy,ixx,iyy,ixxx,iyyy:CARDINAL;
BEGIN
  x:=ABS(x);
  y:=ABS(y);
  ix:=TRUNC(x);
  fx:=x-LFLOAT(ix);
  iy:=TRUNC(y);
  fy:=y-LFLOAT(iy);

  IF fx<0.5 THEN IF ix>0 THEN DEC(ix) END;
  ELSE fx:=fx-1.0 END;
  IF fy<0.5 THEN DEC(iy) ELSE fy:=fy-1.0 END;

  ix:=ix MOD NOISESIZE;
  iy:=iy MOD NOISESIZE;
  ixx:=(ix+1) MOD NOISESIZE;
  iyy:=(iy+1) MOD NOISESIZE;
  ixxx:=(ix+2) MOD NOISESIZE;
  iyyy:=(iy+2) MOD NOISESIZE;

  RETURN qint(
      qint(noisemap[ix  ][iy  ], noisemap[ixx ][iy  ], noisemap[ixxx][iy  ], fx),
      qint(noisemap[ix  ][iyy ], noisemap[ixx ][iyy ], noisemap[ixxx][iyy ], fx),
      qint(noisemap[ix  ][iyyy], noisemap[ixx ][iyyy], noisemap[ixxx][iyyy], fx),fy);
END qrand;

PROCEDURE lrand(x:REAL):REAL;                (* linear interpolated random vector *)
VAR fx:REAL;
    ix:CARDINAL;
BEGIN
  x:=ABS(x);
  ix:=TRUNC(x);
  fx:=x-FLOAT(ix);
  RETURN noisevec[ix MOD (HIGH(noisevec)+1)]*(1.0-fx) + noisevec[(ix+1) MOD (HIGH(noisevec)+1)]*fx;
END lrand;

PROCEDURE dlrand(x:REAL):REAL;                (* derivative linear interpolated random vector *)
VAR ix:CARDINAL;
BEGIN
  x:=ABS(x);
  ix:=TRUNC(x);
  RETURN noisevec[(ix+1) MOD (HIGH(noisevec)+1)] - noisevec[ix MOD (HIGH(noisevec)+1)];
END dlrand;



PROCEDURE rotvecl3(VAR v:LVEC; r:ROTVEC);
VAR h:REAL;
BEGIN
  h  :=v.x*r.cx + v.y*r.sx;
  v.y:=v.y*r.cx - v.x*r.sx;
  v.x:=h  *r.cy + v.z*r.sy;
  v.z:=v.z*r.cy - h  *r.sy;
END rotvecl3;

PROCEDURE rotvec3(VAR v:NORMVEC; r:ROTVEC);
VAR h:REAL;
BEGIN
  h  :=v.x*r.cx + v.y*r.sx;
  v.y:=v.y*r.cx - v.x*r.sx;
  v.x:=h  *r.cy + v.z*r.sy;
  v.z:=v.z*r.cy - h  *r.sy;
END rotvec3;

PROCEDURE rotveclight(VAR xsun:REAL; v:NORMVEC; r:ROTVEC);
BEGIN
  xsun:=(v.x*r.cx + v.y*r.sx)*r.cy + v.z*r.sy;
END rotveclight;


PROCEDURE treepos(pos:POSITIONL; odd:BOOLEAN; VAR th:TREEHINT);
VAR la, lo, sh:LONGREAL;
    dx, dy:REAL;
    onx, ony:INTEGER;
    attr:CARD8;
    stampos:POSITIONL;
BEGIN
  WITH th DO
    sx:=TREERAND*lrand(pos.long*treescramblong);
    sy:=TREERAND*lrand(pos.lat *treescramblat);
    la:=pos.lat *treerasterlat  + sx;
    lo:=pos.long*treerasterlong + sy;

--dxn:=RANDN*dlrand(pos.long*TREESCAMBLONG);
--dyn:=RANDN*dlrand(pos.lat *TREESCAMBLAT);

    IF odd THEN la:=la+0.5; lo:=lo+0.5 END;
--    onx:=VAL(INTEGER, floor(lo+0.25))*2;
--    ony:=VAL(INTEGER, floor(la+0.25))*2;

    hity:=(fracint(la+0.25, ony)-0.5)*2.0;             (* radius from center tree pos / rastersize *)
    hitx:=(fracint(lo+0.25, onx)-0.5)*2.0;
    onx:=onx*2; ony:=ony*2;

    sh:=0.0;
    IF odd THEN onx:=onx+1; ony:=ony+1; sh:=0.5 END; 
    IF (nx<>onx) OR (ny<>ony) THEN                     (* another tree *)

      nx:=onx;
      ny:=ony;
      stampos.lat :=(VAL(LONGREAL, ny)*0.5-sh)/treerasterlat;
      stampos.long:=(VAL(LONGREAL, nx)*0.5-sh)/treerasterlong;

      dx:=TREERAND*lrand(stampos.long*treescramblong);
      dy:=TREERAND*lrand(stampos.lat *treescramblat );
      stampos.lat :=stampos.lat -dx/treerasterlat;
      stampos.long:=stampos.long-dy/treerasterlong;

      normvec(stampos, ADR(meta), tslant, rooth, nvec, water, wood, urban);
--IF debugxi=TESTX THEN WrFixed(stampos.lat/RAD, 6, 15); WrFixed(stampos.long/RAD, 6, 15); WrFixed(rooth, 2, 9); WrStrLn(" rooth ") END; 
      stamtoeye:=distnear(teye, stampos);
    END;
  END; 
END treepos;


PROCEDURE treesizealt(alt:REAL):REAL;
BEGIN RETURN treesize*one((zerotreealt-alt)*treeshrink) (**one(dist*(1.0/SHOWTREEDIST))*) END treesizealt;

PROCEDURE wooddensity(attw:CARD8; alt, slant:REAL):REAL;
--BEGIN RETURN FLOAT(attw)*0.01*one((2100.0-alt)*0.0025)*one((slant-0.4)*botanicslant) END wooddensity;
--BEGIN RETURN (* one(FLOAT(attw)*0.02-1.0)* *) min(one((2100.0-alt)*0.0025),one((slant-0.4)*botanicslant)) END wooddensity;
--BEGIN RETURN min(sqr(one(treesizealt(alt)*0.033)),one((slant-0.4)*botanicslant)) END wooddensity;
BEGIN RETURN one((slant-0.4)*botanicslant) END wooddensity;




PROCEDURE treeform(pos:POSITIONL; eyeh:REAL; odd:BOOLEAN; VAR stacked:BOOLEAN; VAR th:TREEHINT):CARD8;
VAR roundness, thickness, noise:REAL;
    inois:CARDINAL;
    stam:BOOLEAN;
BEGIN
  treepos(pos, odd, th);
  WITH th DO
--IF (nx<>294325) OR (ny<>1519253) THEN (*debugg:=500.0;*) RETURN 0 END;
    IF stamtoeye<SHOWTREEDIST THEN RETURN 0 END;

    inois:=CRC[ny MOD 256 + nx MOD 256*256];
    noise:=FLOAT(inois)*(1.0/65536.0);
--WrInt(wood,3); WrStrLn("");
    IF (wood<50 +CRC[(nx*11+ny*4) MOD 65536] DIV 2048) OR (wooddensity(wood, rooth, tslant)-TREENOTCH<noise) THEN RETURN 0 END;  (* remove random trees *)
--    IF (wood<50) THEN RETURN 0 END;  (* remove random trees *)
--WrInt(nx, 10); WrInt(ny, 10); WrStrLn("nxny");

    regionpeak:=treesizealt(rooth)*one(0.1+plus(tslant-0.4)*botanicslant); (* treesize on altitude and slant *)
--    areapeak:=regionpeak*one(FLOAT(wood)*0.01+0.5);
areapeak:=regionpeak;
    noise:=FLOAT(CRC[(nx*17+ny*3) MOD 65536]);
    peakh:=areapeak*(1.0-noise*(0.4/65536.0));  (* random tree size *)
    IF peakh<MINTREESIZE THEN RETURN 0 END;

--    thickness:=1.0 (*sqr(20.0/peakh)*);
--    thickness:=regionpeak/peakh;
--    hitx:=(pos.long-stampos.long)*treerasterlat*2;
--    hity:=(pos.lat -stampos.lat )*treerasterlong*2;
    hitradius:=hitx*hitx + hity*hity;                         (* radius^2 normalized to treeraster *)
    thickness:=(treerastersize*2.0)/peakh;
    IF thickness<1.0 THEN thickness:=1.0 END;
    thickness:=hitradius*thickness;
    IF thickness<1.0 THEN hith:=peakh*(1.0-pow(thickness, TREEFORM));
    ELSE hith:=0.0 END;
    og:=eyeh-rooth;
    stam:=og<thickness*treerastersize*0.6 + peakh*(noise+20000.0)*(STAMLEN/(20000+65536.0));
    isstam:=stam & (thickness<=STAMTHICKNESS);
    IF isstam & (og<0.0) THEN RETURN 0 END;

    IF NOT isstam & (stamtoeye<showstamdist) & stam THEN       (* under tree but not stam *)
      stacked:=TRUE;
    RETURN 1 END;                                              (* transparent *)

    RETURN 2                                                   (* hit tree *)
  END;
END treeform;



(*
PROCEDURE treevec(pos:POSITIONL; VAR nvec:NORMVEC);
VAR x,y,r:REAL;
BEGIN
  treepos(pos, x, y);
  r:=osi.sqrt(x*x + y*y);
  IF r<>0 THEN r:=0.5/r ELSE r:=0.0 END;  
  nvec.z:=0.85;
  nvec.x:=-x*r;
  nvec.y:=-y*r;
END treevec;	
*)


PROCEDURE raytrace(x0,y0,z0, dx,dy,dz:LONGREAL; maxdist, altrise, refrac, treesize:REAL;
                   mirrord:BOOLEAN; VAR dist:LONGREAL; VAR subpix:REAL; VAR pos:POSITIONL;
                   VAR treehints:TREEHINTS; VAR ddstack:DDSTACK; VAR istree, odd:BOOLEAN;
                   shad:BOOLEAN);
CONST LINSTEP=4000.0;

VAR sp, spo, lastsp, step, resol, ypixalt, tol, tolw, hm, hmm, treestep:REAL;
    iter:CARDINAL;
    attr, od:CARD8;
    meta:METAINFO;
    alt, h, dlat,dlong,dalt, lat0, long0, alt0, lat1, long1, alt1, dist0, dist1, ddist:LONGREAL;
    stk:BOOLEAN;
BEGIN
  step:=0.0;
  IF dist=0.0 THEN lastsp:=0.0 ELSE lastsp:=MAX(REAL) END;

  alt:=MAX(REAL);
  subpix:=0.0;
  spo:=0.0;
--WrStrLn("---+ ");
--WrFixed(dist, 1,5); WrStr("=dist ");
  meta.withslant:=shad;
  meta.aliasattr:=TRUE;                                                (* watch for wood *)
  iter:=5;
--  oldh:=MAX(REAL);
  dist0:=dist;
  dist1:=dist;
  wgs84rl(x0 + dx*dist1, y0 + dy*dist1, z0 + dz*dist1, lat1, long1, alt1);
  alt1:=alt1 - dist1*dist1*refrac;

  IF (ddstack.startdist=dist) & (alt1>ddstack.stoppalt) THEN ddstack.startdist:=0.0 END;

  LOOP
    treehints[FALSE].peakh:=0.0;treehints[TRUE].peakh:=0.0;
    IF dist>=dist1 THEN                                                (* end of linearized peace *)
      lat0:=lat1;
      long0:=long1;
      alt0:=alt1;
      dist0:=dist1;
      dist1:=dist1+LINSTEP;
      wgs84rl(x0 + dx*dist1, y0 + dy*dist1, z0 + dz*dist1, lat1, long1, alt1);
      alt1:=alt1 - dist1*dist1*refrac;
      dalt :=(alt1 -alt0 )*(1.0/LINSTEP);
      dlong:=(long1-long0)*(1.0/LINSTEP);
      dlat :=(lat1 -lat0 )*(1.0/LINSTEP);
    END;
    ddist:=dist-dist0;                                                 (* go ahead linearized *)
    pos.lat :=lat0 +dlat *ddist;
    pos.long:=long0+dlong*ddist;
    alt     :=alt0 +dalt *ddist;
    ypixalt:=ABS(altrise)*dist;                                        (* pixel size *)
    tol:=ypixalt*0.25;
    IF tol<0.1 THEN tol:=0.1 END;
    h:=getsrtmlong(pos.lat, pos.long, TRUNC(ABS(ypixalt)*0.5), ypixalt<5.0, resol, attr, ADR(meta));
    IF defaultocean & ((h>NOALT) OR (h=0.0) & (resol>=900.0)) THEN h:=0.0 END; (* no srtm or srtm30=0 is water *)

    IF h<10000.0 THEN 
      sp:=alt-h;
--IF altrise<0.0 THEN WrFixed(h, 1,7); WrFixed(pos.lat/RAD, 6,12); WrFixed(pos.long/RAD, 6,12); END;
--IF shad & (debugxi=10) THEN WrFixed(alt, 2, 8); WrStr(","); WrFixed(h, 2, 0); END;
      odd:=FALSE;
      istree:=FALSE;
      IF (treesize>0.0) & (sp<=treesize) & NOT alticolour & (meta.attrweights[ATTRWOOD]>0.0) (*& (dist>SHOWTREEDIST)*) THEN
        IF treehints[FALSE].treepass THEN                                           (* use constant wood higth for poi *)
          IF (meta.attrweights[ATTRWOOD]>0.5) & (sp<treesizealt(h))
          & (sp<treehints[FALSE].hith) & NOT mirrord THEN
            treehints[FALSE].hith:=sp;
            treehints[FALSE].passdist:=dist;
          END; 
        ELSE 

          stk:=FALSE;
          odd:=FALSE;
          od:=treeform(pos, alt, FALSE, stk, treehints[FALSE]);
          IF (od=2) & (treehints[FALSE].hith+treehints[FALSE].rooth<alt) THEN od:=0 END;
          IF od<2 THEN
            od:=treeform(pos, alt, TRUE, stk, treehints[TRUE]);

            IF (od=2) & (treehints[TRUE].hith+treehints[TRUE].rooth<alt) THEN od:=0 END;
            IF od>0 THEN odd:=TRUE END;
          END;
          IF stk & NOT mirrord & (ddstack.startdist=0.0) THEN
            ddstack.startdist:=dist;
            ddstack.stoppalt:=treehints[odd].peakh+treehints[odd].rooth;
          END;
          IF (od=2) & (treehints[odd].hith>0.0) THEN
            istree:=TRUE;
            sp:=alt-(treehints[odd].rooth+treehints[odd].hith);      (* add treeform to srtm *)
          END;

        END;
        tolw:=tol*0.25;

      ELSE tolw:=tol END;

      IF shad THEN
--IF (debugxi=10) THEN WrFixed(sp, 2, 8) END;
        IF dist>0.0 THEN
          hmm:=1.0+ABS(meta.slantx)+ABS(meta.slanty);          (* aproximativ distance correction to slants *)
          IF istree & (hmm<10.0) THEN hmm:=10.0 END;
          hm:=hmm*treehints[FALSE].aliasm;                         (* pixl size on shadow point *)
          sp:=sp+hm*0.5;
          IF istree & (hm<5.0) THEN hm:=5.0 END;               (* soft tree shadow *)
          hmm:=dist*(SUNSIZE*2.0);
          IF hmm>hm THEN hm:=hmm END;                          (* sun diameter more important *)
          IF hm<2.0 THEN sp:=sp+0.1*irand(pos.lat*3000000.0, pos.long*coslat*3000000.0) END;
          hmm:=sp/hm;
          IF hmm<treehints[FALSE].missspace THEN       (* store lowest pass angle over obstakle *)
            treehints[FALSE].missspace:=hmm;
          END;
          sp:=sp+hm;
          IF sp<0.0 THEN EXIT END; 
        END;
        IF sp>0.0 THEN
          step:=sp*0.3+0.1;
          IF (step>1.0) & (treesize>0.0) & (sp<=treesize)
          & ((meta.attrweights[ATTRWOOD]>0.0) OR (dist<80.0)) THEN step:=1.0 END;
        ELSE step:=0.2 END; 

--      IF shad THEN AddShadow(dist, h+sp) END;

      ELSIF sp>tolw THEN                                           (* over ground *) 


--        IF (treesize>0.0) & (sp<=treesize*1.2) & (meta.attrweights[ATTRWOOD]>0.0)
--        THEN step:=ypixalt+0.2 ELSE step:=sp*0.3 END;


--        step:=sp*0.3;
        step:=sp*0.5;
        IF spo>0.0 THEN step:=step*sp/spo END;         (* shorter step on fast dropdown *)
        IF (treesize>0.0) & (sp<=treesize*1.2) THEN
          IF (meta.attrweights[ATTRWOOD]>0.0) & (step>ypixalt+0.2) THEN step:=ypixalt+0.2 END;
        END;
--        IF step>treerastersize*0.5 THEN step:=treerastersize*0.5 END;




        IF subpix=0.0 THEN 
          IF sp<lastsp THEN lastsp:=sp;
          ELSIF ypixalt<>0.0 THEN subpix:=lastsp/ypixalt END;
        END;
      ELSIF sp>-tolw THEN EXIT                                     (* hit earth *)

      ELSE                                                         (* under ground *)
IF mirrord THEN EXIT END;                                   (* quick *) 
        IF spo<=0.0 THEN                                           (* again under ground *)
          step:=sp*0.5;
--IF mirrordist>0.0 THEN step:=step*0.1 END;

--WrStr(" back-back ");
--EXIT; 
        ELSE   
          DEC(iter);
          IF iter=0 THEN EXIT END;

          step:=sp/(spo-sp)*step;           (* first time under ground jump back with newton *)
        END;
--WrFixed(step, 1,6); WrFixed(spo-sp, 1,7); WrInt(iter, 3); WrStrLn(" newton ");

      END;
--      oldh:=h; 
    ELSIF ABS(step)<tolw THEN                                    (* jump over hole in erth *)
      DEC(iter);
      IF iter=0 THEN EXIT END;

      sp:=0.0; 
      IF step<30.0 THEN step:=30.0 END;
    END;
    spo:=sp;
    dist:=dist + step;
    IF (dist<-10.0) OR (dist>maxdist) OR (alt>maxmountain) & (dalt>0.0) THEN dist:=maxdist+100.0; EXIT END;   (* we are in dust or heaven *)
  END;

  IF subpix>1.0 THEN subpix:=1.0 END;
END raytrace;
(* antialiasing: search for peaks and set highest pixel smooth *)



PROCEDURE rotvector(VAR a,b:LONGREAL; cw,sw:LONGREAL);
VAR h:LONGREAL;
BEGIN
  h:=a*cw+b*sw;
  b:=b*cw-a*sw;
  a:=h
END rotvector;



PROCEDURE angtoxyz(wx,wy,slat,clat,slong,clong:LONGREAL; VAR x,y,z:LONGREAL);
VAR cwy:LONGREAL;
BEGIN
  cwy:=cos(wy);
  x:=-sin(wy);
  y:=sin(wx)*cwy;
  z:=cos(wx)*cwy;
  rotvector(z,x, clat,slat);
  rotvector(y,x, clong,slong);
  x:=-x;
END angtoxyz;


PROCEDURE posdirtoxyz(pos-:POSITIONL; az,el:LONGREAL; VAR xn,yn,zn:LONGREAL);
BEGIN
  angtoxyz(az,el, sin(-pos.lat),cos(-pos.lat), sin(-pos.long),cos(-pos.long), xn,yn,zn);
END posdirtoxyz;

PROCEDURE rotateview(VAR wx, wy:LONGREAL; xim, yim, incl:LONGREAL);
VAR z:LONGREAL;
BEGIN
  z:=1.0;
  rotvector(yim,z,cos(incl),sin(incl));
  wx:=atang2l(z, xim);
  wy:=atang2l(sqrt(xim*xim + z*z), yim);
END rotateview; 


PROCEDURE seesun(frompos-:POSITIONL; alt:LONGREAL; sunaz, sunele, maxdist, refrac, pixperm:REAL):REAL;
VAR dd, xn,yn,zn, x0,y0,z0:LONGREAL;
    pos:POSITIONL;
    th:TREEHINTS;
    space:REAL;
    istree, odd:BOOLEAN;
    ddstack:DDSTACK;
BEGIN
  posdirtoxyz(frompos, sunaz, sunele, xn,yn,zn);
  wgs84sl(frompos.lat, frompos.long, alt, x0,y0,z0);
  FILL(ADR(th), 0C, SIZE(th));
  odd:=FALSE;
  dd:=0.0;
  th[FALSE].missspace:=MAX(REAL);
  th[FALSE].aliasm:=1.0/pixperm;
  ddstack.startdist:=0.0;
  raytrace(x0,y0,z0, xn,yn,zn, maxdist, 0.0, refrac, treesize, TRUE, dd, space, pos, th, ddstack,
           istree, odd, TRUE);

  IF dd>maxdist THEN RETURN one(th[FALSE].missspace+0.5) END; (* sight to sun *)

  RETURN 0.0;
END seesun;


PROCEDURE sundisk(x, wy, sunazi, sunele:REAL; dist, sunarea, maxdist:REAL; lenseffect:BOOLEAN;
                  suncol-:COL; refrac:REAL; VAR c:COL);  (* dust and all else in heaven *)

  PROCEDURE inv(x:REAL):REAL;
  BEGIN
    x:=ABS(x);
    IF x=0.0 THEN RETURN MAX(REAL) END;
    RETURN 1.0/(x)
  END inv;

CONST MAXSUNDUST=3600.0;  (* 4000 *)
      LENSEFFECT=40.0;
      HEAVENDUST=1.0;     (* 2 *)
      RAD90=1.0/(90*RAD);
       
VAR y, sd, s, sy, ed, as, ash, lens, bl:REAL;
    hc,nh:COL;

BEGIN 
  y:=radmod(wy-sunele);
  sd:=x*x+y*y;
  bl:=maxdist;
  IF dist<>0.0 THEN bl:=dist END;
  bl:=bl*(RAD90/15000);                                 (* lesser blue sky on shorter sight *)
  IF bl>RAD90 THEN bl:=RAD90 END;
  nh.r:=1.0;
  nh.g:=1.0;
  nh.b:=1.0;
  hc.r:=0.7;
  hc.g:=1.0;
  hc.b:=1.5;
  colblend(bl, hc, nh, hc);                             (* near water dust to distant oxygene blue *)
  colblend(ABS(x)*0.3, suncol, hc, hc);                 (* azimuth sun colour in dust *)
 
  IF dist<>0.0 THEN                                     (* earth *)
    colblend((0.5-sqrt(sd)*3.0)*sunarea, hc, suncol, hc); (* sun angle dust colour on earth *)
  END;


  s:=600.0;
  IF sunele<0.0 THEN                                    (* night *)
    s:=s+sunele*3500.0;
    IF s<20.0 THEN s:=20.0 END;
  END;

  hc.r:=s*hc.r;
  hc.g:=s*hc.g;
  hc.b:=s*hc.b;
        
  IF dist=0.0 THEN                                      (* heaven *)
    IF wy<0.0 THEN wy:=0.0 END;
    ed:=1.0-wy*bl;
    s:=HEAVENDUST+x*x+ABS(sunele);                             
    IF s>6.0 THEN s:=6.0 END;                           (* dust elevation *)
    IF ed<=0.0 THEN ed:=0.0 ELSE ed:=pow(ed, s) END;    (* ray elevation effect *)

    IF sunele<0.0 THEN                                  (* night heaven blue *)
      s:=1.0+sunele*10.0;
      IF s<0.02 THEN s:=0.02 END;
      nh.r:=heavencol.r*s;
      nh.g:=heavencol.g*s;
      nh.b:=heavencol.b*s;
    ELSE nh:=heavencol END;
    colblend(ed, nh, hc, c);                            (* blue sky down to dust *) 

  ELSE                                                  (* earth *)
    s:=1-dist; 
    c.r:=c.r*s;
    c.g:=c.g*s;
    c.b:=c.b*s;
    colblend(dist, c, hc, c);
  END;

  IF lenseffect THEN                                    (* direct sun *)
    IF sd<=SUNSIZE*SUNSIZE THEN                         (* view into sun *)
      lens:=SUNLUM;
    ELSE
      lens:=0.5+inv(ABS(sin(atang2(x,y)*7.0)));
      IF lens>1.5 THEN lens:=1.5 END;
      lens:=lens/(sd*250.0+0.1);
      lens:=lens*sunarea*LENSEFFECT;
    END;

    c.r:=c.r+suncol.r*lens;
    c.g:=c.g+suncol.g*lens;
    c.b:=c.b+suncol.b*lens;
  END;

END sundisk;

PROCEDURE waterglitter(x, wy, sunazi, sunele:LONGREAL; pixsize, maxdist:REAL;
                  suncol-:COL; pos-:POSITIONL; refrac, alt, dropsun:REAL; VAR c:COL);  (* small waves mirrors sun on water *)
CONST ZOOM=4000.0;
VAR y, lens:LONGREAL;
BEGIN 
  y:=radmod(wy-sunele);
  lens:=(2.0-sunele)*0.12/(x*x+y*y*sunele*(1.0/PI)+0.0005)*(1.0+sqr(qrand((x+1.0)*ZOOM+FLOAT(millisec)*0.0003, (y+1.0)*ZOOM)*(1.0-one(pixsize*0.2))));
  IF lens>0.1 THEN
    IF NOT dropshadows THEN
      dropsun:=seesun(pos, alt+WATERMIRRORGAP, sunazi, sunele+SUNSIZE, maxdist, refrac, 1.0);
    END;
    IF dropsun>0.0 THEN
      lens:=lens*dropsun*waterglance; 
      c.r:=c.r+suncol.r*lens;
      c.g:=c.g+suncol.g*lens;
      c.b:=c.b+suncol.b*lens;
    END; 
  END;
END waterglitter;



PROCEDURE sunvisability(sunaz, sunel,slat,clat,slong,clong, x0,y0,z0:LONGREAL; maxdist, refrac:REAL):REAL;
CONST STEPS=8;
VAR dd, wx, wy, xn,yn,zn:LONGREAL;
    space, sum,r:REAL;
    pos:POSITIONL;
    i:CARDINAL;
    th:TREEHINTS;
    istree, odd:BOOLEAN;
    ddstack:DDSTACK;
BEGIN
  r:=0.0;
  sum:=0.0;
  FILL(ADR(th), 0C, SIZE(th));
  odd:=FALSE;
  FOR i:=0 TO STEPS-1 DO
    wx:=sunaz+(SUNSIZE*0.5)*osi.sin(r);
    wy:=sunel+(SUNSIZE*0.5)*osi.cos(r);
    angtoxyz(wx,wy,slat,clat,slong,clong, xn,yn,zn);
    dd:=0.0;
    ddstack.startdist:=0.0;
    raytrace(x0,y0,z0, xn,yn,zn, maxdist, 0.0, refrac, treesize, FALSE, dd, space, pos, th, ddstack,
    istree, odd, FALSE);
--WrFixed(sum, 2, 7);
    IF dd>maxdist THEN sum:=sum+1.0/STEPS END;
    r:=r+(PI*2.0/STEPS);
  END;
--WrStrLn(" sums");
  RETURN sum
END sunvisability;



PROCEDURE snowcol(VAR c:COL; sun-:COL; xsun, sunlum:REAL);
VAR x:REAL;
BEGIN
  IF xsun<-0.3 THEN xsun:=-0.3 END;
  x:=(xsun+2.0)*sunlum;
--  x:=xsun+2.0;
  c.r:=cSNOWSHADOWr*x;
  c.g:=cSNOWSHADOWg*x;
  c.b:=cSNOWSHADOWb*x;
  IF xsun>0.0 THEN                      (* surface sees sun *)
    c.r:=c.r + cSNOWr*xsun*sun.r;
    c.g:=c.g + cSNOWg*xsun*sun.g;
    c.b:=c.b + cSNOWb*xsun*sun.b;
  END; 
END snowcol;

PROCEDURE botaniccol(VAR c:COL; suncol-:COL; suncos, diffuslum:REAL);
VAR x:REAL;
BEGIN
--  x:=(xsun+2.5)*(sun.r+sun.g+sun.b)*0.33333333;
  IF suncos<-0.3 THEN suncos:=-0.3 END;
  x:=(suncos+3.5)*diffuslum;
  c.r:=0.05*x;
  c.g:=0.055*x;
  c.b:=0.09*x;
  IF suncos>0.0 THEN 
    c.r:=c.r + suncos*suncol.r;
    c.g:=c.g + suncos*suncol.g;
    c.b:=c.b + suncos*suncol.b;
  END; 
END botaniccol;

PROCEDURE treecol(pos:POSITIONL; th-:TREEHINT; sunv-:ROTVEC;
                  pixsize, diffuslum, dense, snowdense, eyesunaz, xsun, eyeang:REAL; hitwood:CARD8;
                  suncol-:COL; VAR c:COL);
VAR h,sh,sun, rnd, trestruc, ang, geosun, pixs, pixst, pixshh, pixsun, hxx, hxy:REAL;
    v:NORMVEC;
    ri, gi, bi:REAL;
    s, ss:COL;
BEGIN
  pixst:=pixsize*0.5;  (* structure amplitude *)
  IF pixst>1.0 THEN pixst:=1.0 END;
  pixs:=one(pixsize*0.5);
  pixsun:=pixsize*1.5;
  IF pixsun>1.0 THEN pixsun:=1.0 END;
  pixshh:=one(pixsize*2.0);
  IF dropshadows THEN geosun:=diffuslum ELSE geosun:=one(diffuslum*20.0) END;

  ang:=PI*1.5+atang2(th.hitx - TREERAND*dlrand(pos.lat *treescramblat)*th.hity,
                 th.hity - TREERAND*dlrand(pos.long*treescramblong)*th.hitx);
  rnd:=irand(FLOAT((th.nx*7+th.ny*6) MOD 512)+ang*5.0, (th.og+30.0)*6.0);   (* tree structure *)
  trestruc:=(1.0+rnd*(1.0/2.0))*pixst + (1.0-pixst);                  (* flatten struct on distance *)

  IF th.hitradius<1.0 THEN
    v.x:=-th.hitx; v.y:=-th.hity; v.z:=1.0-osi.sqrt(th.hitradius)*treesize/th.peakh;
    rotveclight(sun, v, sunv);                                   (* sun angle effekt *)
  ELSE sun:=-1.0 END;


  IF xsun>0.0 THEN
    xsun:=xsun*(0.75+ABS(eyesunaz)*(0.5/PI)); (* distant in suncos median wood brightness *)
  END;
  sun:=pixsun*sun + (1.0-pixsun)*xsun;        (* lower individual sun angle on dist *)

  IF dropshadows THEN sh:=1.0 ELSE sh:=(sun*0.15+0.85) END;                        (* sun angle effekt *)

  IF th.isstam THEN                                              (* stam *)
    snowdense:=0.0;
    h:=irand(FLOAT(th.nx*7+th.ny*5), FLOAT(th.nx*4+th.ny*3))*cSTAMVAR; 
    sh:=sh*trestruc;
    ri:=sh*cSTAMr;
    gi:=sh*(cSTAMg+h);
    bi:=sh*(cSTAMb+h);
  ELSE
    snowdense:=one(snowdense*1.6-pixs*ABS(2.0-ABS(qrand(FLOAT((th.nx*13+th.ny*6) MOD 512)+ang*2.0, th.og*4.0))));
    snowdense:=snowdense*pixs + snowdense*qrand(FLOAT(th.nx), FLOAT(th.ny))*(0.1-pixs)*pixshh + snowdense*(0.35-pixshh*0.35);

    IF dropshadows THEN
      h:=th.og/th.peakh;                                  (* more light on higth *)
      IF h<0.3 THEN h:=0.3 END;      
      h:=(h+0.2);
    ELSE
      h:=sh*(1.0-one((0.8-th.hith/th.areapeak)*2.5)*1.2*one(FLOAT(hitwood)*0.02-1.0)*pixs);  (* use relative hight dark in wood *) 
      h:=h+pixs*one((th.hith/th.regionpeak-0.8)*4.0)*0.6;            (* illuminate high peaks *) 
    END;
--    h:=(1.0+rnd*(1.0/4.0))*pixst*h (* + ((1.0-eyeang*0.4)-pixst) *);
    h:=h*trestruc;
    IF h<0.1 THEN h:=0.1 END;
    ri:=cTREEr*h;
    gi:=cTREEg*h;
    bi:=cTREEb*h;
  END;

  ri:=ri*(1.0+irand(FLOAT(th.nx)*0.2, FLOAT(th.ny)*0.2)*0.04);   (* colour variance of trees *)
  gi:=gi*(1.0+irand(FLOAT(th.ny)*0.2, FLOAT(th.nx)*0.2)*0.08);

  botaniccol(c, suncol, sun, diffuslum);
  c.r:=c.r*ri;
  c.g:=c.g*gi;
  c.b:=c.b*bi;


  IF snowdense>0.0 THEN
    snowcol(s,  suncol, sun, geosun);
    colblend(snowdense, c, s, c);             (* mix green with snow *)
  END;
END treecol;


PROCEDURE citycol(attu:CARD8; pos-:POSITIONL; pixs, viewang:REAL; VAR c, lamp:COL);
CONST PARCEL=210000.0;
VAR street,hr,hw,gr, gw, gravel,ax,ay,lx,ly,l,vis,sh, med, str:REAL;
    ry, rx:INTEGER;
    hc, hh:COL;
    i:CARDINAL;
BEGIN
--  viewang:=(ABS(viewang)+0.5)*pixs*10.0;
  viewang:=15.0*pixs; 
  lamp.r:=0.0;lamp.g:=0.0;lamp.b:=0.0;
  IF (attu>0) & (urbanoff<>0.0) THEN
    lx:=fracint(pos.lat*PARCEL, rx);
    ly:=fracint(pos.long*coslat*PARCEL, ry);
    ax:=lx;
    ay:=ly;
    IF ax>0.5 THEN ax:=1.0-ax END;
    IF ay>0.5 THEN ay:=1.0-ay END;

    street:=min(ax, ay);
    hw:=10.0*irand(LFLOAT(rx),LFLOAT(ry));
    hr:=0.125*frac(hw);
    vis:=FLOAT(attu)*0.01*one(4.0-pixs/urbanoff*4.0);

    str:=one((0.05-street)*viewang+0.5);
    med:=1.0-str;
    hc.r:=cCITYMEADOWr*med + cCITYSTREETr*str;
    hc.g:=cCITYMEADOWg*med + cCITYSTREETg*str;
    hc.b:=cCITYMEADOWb*med + cCITYSTREETb*str;

    gravel:=irand(pos.lat*30000000.0, pos.long*coslat*30000000.0);
    gr:=one(hw-8.0);
    gw:=one(hw+3.0);
    gravel:=1.0-(1.0-one(sqrt(ABS(gravel))))*0.8*(gw-gr)*one(pixs*0.1);
    hh.r:=cCITYCOL*(0.8+gw)*gravel;
    hh.g:=cCITYCOL*(1.0+gw-gr)*gravel;
    hh.b:=cCITYCOL*(1.0+gw-gr)*gravel;

    colblend((street-(0.15+hr))*viewang+0.5, hc, hh, hc);


(* 
    IF street<0.07 THEN hc.r:=cCITYMEADOWr; hc.g:=cCITYMEADOWg; hc.b:=cCITYMEADOWb;         (* meadow *)
    ELSIF street<0.23+hr THEN hc.r:=cCITYSTREETr; hc.g:=cCITYSTREETg; hc.b:=cCITYSTREETb;   (* street *)
    ELSIF street<0.3+hr  THEN hc.r:=cCITYMEADOWr; hc.g:=cCITYMEADOWg; hc.b:=cCITYMEADOWb;   (* meadow *)
    ELSE                                                   (* house *) 
      gravel:=irand(pos.lat*30000000.0, pos.long*coslat*30000000.0);
      hr:=one(hw-8.0);
      hw:=one(hw+3.0);
      gravel:=1.0-(1.0-one(sqrt(ABS(gravel))))*0.8*(hw-hr)*one(pixs*0.1);
      hc.r:=cCITYCOL*(0.8+hw)*gravel;
      hc.g:=cCITYCOL*(1.0+hw-hr)*gravel;
      hc.b:=cCITYCOL*(1.0+hw-hr)*gravel;
    END;
*)
--lights
    i:=CRC[VAL(CARDINAL,ry) MOD 256 + VAL(CARDINAL,rx) MOD 256*256];
    sh:=FLOAT(1 + i DIV 32 MOD 8); 
    CASE i MOD 32 OF
      0    :lamp.r:=cLAMPTr; lamp.g:=cLAMPTg; lamp.b:=cLAMPTb; sh:=1.0;
     |1..7 :lamp.r:=cLAMPHr; lamp.g:=cLAMPHg; lamp.b:=cLAMPHb;
     |8..15:lamp.r:=cLAMPWr; lamp.g:=cLAMPWg; lamp.b:=cLAMPWb;
    ELSE END;
    l:=vis*sh/(sqr(lx-FLOAT(i DIV 1007 MOD 4 + 1)*0.2)+sqr(ly-FLOAT(i DIV 771 MOD 4 + 1)*0.2) + 0.002*sh);
    lamp.r:=lamp.r*l; lamp.g:=lamp.g*l; lamp.b:=lamp.b*l;
---lights

    colblend(vis, c, hc, c);
  END;
END citycol;


PROCEDURE citylight(VAR c:COL; col, boc, citylamp:COL; attu:CARD8);
VAR dust:REAL;
BEGIN
  dust:=FLOAT(attu);
  c.r:=col.r*boc.r + citylamp.r + dust*cCITYDUSTr;
  c.g:=col.g*boc.g + citylamp.g + dust*cCITYDUSTg;
  c.b:=col.b*boc.b + citylamp.b + dust*cCITYDUSTb;
END citylight;

PROCEDURE snownoise(pos-:POSITIONL; sunang, pixs:REAL):REAL;
VAR h,s,f:REAL;
BEGIN
  f:=qrand(pos.long*40000000.0, pos.lat*40000000.0);
  f:=f*f*f*0.2;
  IF f>1.0 THEN f:=1.0 ELSIF f<-1.0 THEN f:=-1.0 END;
  f:=f*one(pixs*0.3);
  h:=qrand(pos.long*6000000.0*coslat, pos.lat*6000000.0);  (* near *)
  h:=h*h*h*0.1;
  IF h>1.0 THEN h:=1.0 ELSIF h<-1.0 THEN h:=-1.0 END;
  h:=f+h*one(pixs*1.5)*plus(1.0-pixs*0.3);
  s:=qrand(pos.long*520000.0*coslat, pos.lat*520000.0);    (* far *)
--  h:=h+s*0.1;
  h:=sunang-0.1+s*0.02+h/(ABS(sunang)+1.0)*0.1;
  IF h>1.0 THEN h:=1.0 ELSIF h<-1.0 THEN h:=-1.0 END;
  RETURN h
END snownoise;


PROCEDURE meadowcol(pos:POSITIONL; slant, alt, pixs, sunang:REAL; VAR sunvec:REAL; VAR c:COL; coast:REAL);
CONST D=2.0;   (* find best distance feeling *)
VAR hh,s,brown,k,d,f,m,n,sv,wet:REAL;
    hc:COL;
BEGIN
  pos.long:=pos.long*coslat;
  f:=one(pixs*0.2)*qrand(pos.long*(30000000.0*D), pos.lat*(30000000.0*D));
  s:=qrand(pos.long* (8000000.0*D), pos.lat* (8000000.0*D));
  m:=one(pixs*2.0)*qrand(pos.long* (2200000.0*D), pos.lat* (2200000.0*D));
  k:=qrand(pos.long*  (400000.0*D), pos.lat*  (550000.0*D));
  n:=qrand(pos.long*   (70000.0*D), pos.lat*   (70000.0*D));
  s:=s+f;

  sv:=1.0-0.25*f;

--  brown:=one(s*4.0-1.0)*(1.0-ops) + one(m*4.0)*ops;         (* high altitude green/yellow noise *)
  brown:=one(1.0-slant*slant                                                   (* slant brown *)
         + one(alt*(0.00025+ABS(pos.lat*0.0003))+(alt+1000.0-snowalt)*0.001)   (* altitude brown *)
         + one((alt-snowalt+slowaltslant*2.0)/slowaltslant)*0.6);              (* snow disturb *)
--  q:=h*h;
--  q:=q*q;
  d:=(ABS(k)-ABS(m))*0.08*(1.7-one(pixs*2.0));                          (* far variance *)
  sunvec:=s*one(pixs*2.0)*0.8;
(*
  c.r:=cMEADOWr*sv*plus(1.0 + brown*0.2 - d - (h+q)*0.05*ops);
  c.g:=cMEADOWg*sv*plus(1.0 - brown*0.12 + d + (s*s*s*0.01 + s*0.015)*ops);
  c.b:=cMEADOWb*sv*plus(1.0 + brown*0.15 - q*0.05*ops);
*)
  wet:=(coast+50.0)*0.02;                                     (* near water green *)
--wet:=0.0;
  IF wet>1.0 THEN wet:=1.0 END;
  c.r:=cMEADOWr*sv*plus(1.0 + brown*0.2 - d - wet*0.2);
  c.g:=cMEADOWg*sv*plus(1.0 - brown*0.12 + d + n*(1.0-brown)*0.1 + wet*0.15);
  c.b:=cMEADOWb*sv*plus(1.0 + brown*0.15);
  hc.r:=cMEADOWVOIDr;
  hc.g:=cMEADOWVOIDg;
  hc.b:=cMEADOWVOIDb;
--  colblend(one(k*k*10.0)+ABS(s)  (* *one(0.25+pixs*0.25)*), hc, c, c);               (* meadow voids *)
--  colblend(one(k*k*(1.0+20.0*ABS(m)))+ABS(s)  (* *one(0.25+pixs*0.25)*), hc, c, c);               (* meadow voids *)
--  colblend(one(ABS(k)*2.0)+ABS(s), hc, c, c);               (* meadow voids *)
  colblend(one(k*k*10.0+m)+ABS(s)+one(0.75-pixs*0.1), hc, c, c);               (* meadow voids *)
--  IF attwater>0 THEN                                          (* near water sand *)
    hh:=one(1.9-wet);
    hh:=hh*hh*(0.5+sv*0.5);
    hc.r:=cCOASTSANDr*hh;
    hc.g:=cCOASTSANDg*hh;
    hc.b:=cCOASTSANDb*hh;
    colblend((wet-0.85+s*0.06)*25.0, c, hc, c);
--  END;
END meadowcol;

PROCEDURE desert(pos-:POSITIONL; pixs:REAL; VAR c:COL);
VAR n:REAL;
BEGIN
   n:=qrand(pos.long*8000000.0, pos.lat*8000000.0);
   n:=n + one(pixs*0.4) * qrand(pos.long*30000000.0, pos.lat*30000000.0);
   n:=n + one(pixs)*qrand(pos.long*450000.0, pos.lat*450000.0);
   n:=n + qrand(pos.lat*117000.0, pos.long*123000.0);
   n:=0.05+plus(1.0+n*0.05*sandgranularity);

   c.r:=desertcol.r*n;
   c.g:=desertcol.g*n;
   c.b:=desertcol.b*n;
END desert;

PROCEDURE rockcol(pos:POSITIONL; slant, alt, pixs, sunang:REAL; attw:CARD8;
                  istree:BOOLEAN; (*glance-:GLANCESET;*) VAR lastsunv:REAL; VAR c:COL);
VAR nr,ns,n1,n2,na,su,rock,nh:REAL;
    hc:COL;
BEGIN
  pos.long:=pos.long*coslat;
  n1:=qrand(pos.long*9000000.0, pos.lat*9000000.0)*one(pixs);
  nr:=n1+qrand(pos.long*110000.0, pos.lat*11000.0);
  rock:=one(1.1*one(plus((alt+slowaltslant*0.5*FLOAT(ORD(istree))-2200.0)*0.005))
        + plus(1.0-plus(slant-0.4-FLOAT(100-attw)*0.002+nr*0.1)*botanicslant)-0.1);  (* rock / green *)
  rock:=rock+rock*0.3*one(qrand(pos.long*1500000.0, pos.lat*1500000.0)*5.0-2.0);     (* rock / green noise *)

  IF (rock>0.0) THEN                                                 (* else save time *)
    IF desertcol.r<0.0 THEN                                          (* not user defined colour *)
      nh:=qrand(pos.long*8000000.0, pos.lat*8000000.0);
      n2:=nh + one(pixs*0.4) * qrand(pos.long*30000000.0, pos.lat*30000000.0);
      ns:=qrand(pos.long*450000.0, pos.lat*450000.0);
      na:=qrand(alt*3.0+nr*5.0, pos.lat*117000.0+pos.long*123000.0); (* altitude dependent structure *)  
      na:=one(na*10.0)*0.3*ABS(ns)*one(pixs);
      n2:=(n2 + ns*0.5)*one(pixs*2);
      na:=na-(ns*0.3+nh*nh*0.7*one(pixs*0.5))*plus(1.0-pixs*3.0);      (* far rock structure *)
      hc.r:=cROCKr*(na+plus(0.95+one(n1*n1*n1-0.1)*0.15 + n2*0.07));
      hc.g:=cROCKg*(na+one(0.95+n2*0.1));
      hc.b:=cROCKb*(na+one(0.95+n2*0.15));
    ELSE
      desert(pos, pixs, hc);
      IF sandgranularity<0.0 THEN
        rock:=1.0;
--      addglance(pos, glance, ABS(sandgranularity), sandgranularity, sunang, pixs, hc);
      END;
      rock:=rock*1.1;
--addglance(pos, glance, ABS(sandgranularity), sandgranularity, sunang, pixs, hc);
      n2:=0.0;
    END;
    colblend(rock, c, hc, c);
    lastsunv:=rock*n2*0.4 + (1-rock)*lastsunv;                       (* blend meadow/rock sun vector variance *)
  END;
  lastsunv:=sunang+lastsunv*0.05;
  IF desertcol.r<0.0 THEN lastsunv:=sunang*(0.95+0.005/(sqr(sunang)+0.005)) END;
END rockcol;

PROCEDURE snowdensity(pos:POSITIONL; slant, alt, pixs:REAL; th-:TREEHINT; attw:CARD8):REAL;
VAR s,h:REAL;
BEGIN
  h:=0.75-sqr(qrand(pos.long*1000000.0*coslat, pos.lat*1000000.0));
  h:=h*one(pixs)*0.2+qrand(pos.long*10000000.0*coslat, pos.lat*10000000.0)*0.2*one(pixs*0.4);
  h:=h*(0.7-newsnow*0.5);
  IF (treesize>0.0) & (attw>0) & (th.peakh>0.0) THEN
    h:=h-0.25*one((1.0-sqrt(th.hitradius))*2.0);                     (* less snow under tree *)
  END;
  s:=(alt-snowalt)/slowaltslant;
  s:=one((s + slant*slant)*slowaltslant*0.005 + 0.5 + h);
  RETURN s*one((pow(slant,0.5)-0.5)*5.0)                             (* never snow on big slant *)
END snowdensity;

PROCEDURE darkundertree(attw:CARD8; wooddens, snowdense, pixsize:REAL; th-:TREEHINT; VAR c:COL);   (* dark under tree *)
VAR a,p:REAL;
    hc,hh:COL;
BEGIN
(*
  IF (treesize>0.0) & (attw>0) THEN                                  (* tree enabled and tree area *)
    a:=1.3-one(FLOAT(attw)*one(wooddens)*treesize*0.0005);
    IF th.peakh>0.0 THEN a:=min(a, 0.3+0.7*sqrt(th.hitradius)) END;  (* dark disk under tree *)
    hc.r:=c.r*0.15;
    hc.g:=c.g*0.20;
    hc.b:=c.b*0.12;
    hh.r:=c.r*0.4;
    hh.g:=c.g*0.4;
    hh.b:=c.b*0.4;
    colblend(snowdense, hc, hh, hc); (* snow or not under trees *)
    colblend(a, hc, c, c);
  END;
*)
  IF (treesize>0.0) & (attw>0) THEN                                  (* tree enabled and tree area *)
--    a:=one(1.3-one(FLOAT(attw)*wooddens*treesize*0.0005));   (* diffuse light reduction near trees *)
    a:=FLOAT(attw)*0.01;                                       (* diffuse light reduction near trees *)
    IF th.peakh>0.0 THEN                                       (* dark disk under tree *)
      a:=max(a, 1.0-sqrt(th.hitradius));
      c.g:=c.g*(0.8+th.hitradius*0.2);                         (* browner near stam *)
      c.b:=c.b*(0.6+th.hitradius*0.4);                         (* browner near stam *)
    END;
--    IF th.peakh>0.0 THEN a:=max(a, 1.0-th.hitradius) END;      (* dark disk under tree *)
    a:=a*wooddens;
    p:=3.0*pixsize;
    IF p>1.0 THEN p:=1.0 END;
    a:=a*p;
    hc.r:=c.r*0.2;                                                   (* far wood colour *)
    hc.g:=c.g*0.22;
    hc.b:=c.b*0.16;
    hh.r:=c.r*0.4;
    hh.g:=c.g*0.4;
    hh.b:=c.b*0.4;
    colblend(snowdense, hc, hh, hc);                                 (* snow or not under trees *)
    colblend(a, c, hc, c);
--IF debugxi=TESTX THEN WrFixed(wooddens, 2, 6) END;
  END;

END darkundertree;

PROCEDURE fogsun(sun, sight, sunel, alt:REAL):REAL;  (* sun thru fog *)
VAR a:REAL;
BEGIN
  alt:=0.7+alt*(0.3/1500);  IF alt>1.0 THEN alt:=1.0 END;
  RETURN sun*one((ABS(sunel)+0.3)*sight*alt*0.0001);
END fogsun;

PROCEDURE underwatercol(pos:POSITIONL; pixs, wsh, snow, wy, wavx, wavy:REAL; sun-:COL; VAR c:COL);
VAR h:REAL;
    hc:COL;
BEGIN
  IF snow>2500.0 THEN snow:=2500.0 END;
  IF watercol.r<0.0 THEN                            (* automatic water colour *)
    c.r:=cWATERALIASr-snow*0.06;
    c.g:=cWATERALIASg-snow*0.04;                    (* guess temperature and botanic activity *)
    c.b:=cWATERALIASb-snow*0.07;
  ELSE c:=watercol END;

--  h:=1.0-ABS(wy)*0.8;                             (* less mirror on high view see deeper in darker *)
  h:=sin(wy);                                       (* less mirror on high view see deeper in darker *)
--  IF h<0.3 THEN h:=0.3 END;
  c.r:=c.r*(1.0-h*0.3);
  c.g:=c.g*(1.0-h*0.4);
  c.b:=c.b*(1.0-h*0.55);

  IF wsh<1.0 THEN                                   (* see ground *)
    IF millisec>0 THEN                              (* else save time *)
      h:=plus(wsh-0.45)*0.00001; 
      pos.long:=pos.long + wavx*h;                  (* water ground move with waves *) 
      pos.lat :=pos.lat  + wavy*h;
    END; 
    h:=one(pixs)*qrand(pos.long* 8000000.0, pos.lat* 8000000.0);
    h:=h+one(pixs*0.2)*qrand(pos.long*25000000.0, pos.lat*25000000.0);
    h:=one(h*h);
    hc.r:=cUNDERWATERr*(1.0+h*0.4);
    hc.g:=cUNDERWATERg*(1.0+h*0.5);
    hc.b:=cUNDERWATERb;

    h:=plus(1.0+wavx*10.0);                         (* wave lens for light on ground *)  
    hc.r:=hc.r*h;
    hc.g:=hc.g*h;
    hc.b:=hc.b*h;

    colblend((1.0-wsh)*2.0*(0.1+ABS(wy))*watertransparency, c, hc, c); (* remove ground colour on deep water *)
  END;
(*
  h:=0.1+ABS(wy)*2.5;                                (* see into water *)
  IF h>1.0 THEN h:=1.0 END;
  c.r:=c.r*h;
  c.g:=c.g*h;
  c.b:=c.b*h;
*)
  IF surfacedust>0.0 THEN                            (* add white on flat view from dust on surface *)
    h:=0.3+snow*0.0003;
    IF h>1.0 THEN h:=1.0 END;
    h:=(1.0-pow(ABS(wy), 0.05))*h*surfacedust;
    IF h<0.0 THEN h:=0.0 END;
    c.r:=c.r+h;
    c.g:=c.g+h;
    c.b:=c.b+h;
  END;

  c.r:=c.r*sun.r;
  c.g:=c.g*sun.g;
  c.b:=c.b*sun.b;
END underwatercol;

PROCEDURE waveang(pos:POSITIONL; waveamp, wavelenx, waveleny, pixs:REAL; VAR wavx, wavy:REAL);
VAR t,r,u,v:LONGREAL;
    a,x,y:REAL;
    i:CARDINAL;
BEGIN
  pos.long:=pos.long*VAL(LONGREAL, waveleny);
  pos.lat :=pos.lat *VAL(LONGREAL, wavelenx);
  IF millisec>0 THEN                               (* save time if not movie *)
    t:=LFLOAT(millisec)/VAL(LONGREAL,(waveamp+1.0)*500.0) + pos.long*(1000000.0/WAVLEN);     (* speed depending on amplitude *)
    i:=TRUNC(t);
    u:=LFLOAT(i);
    r:=t-u;
    IF ODD(i) THEN r:=1.0-r END;
    u:=LFLOAT(i DIV 2)*10.0;
    v:=LFLOAT((i+1001) DIV 2)*10.0;
    x:=qrand(pos.long+v, pos.lat)*(1.0-r) + qrand(pos.long, pos.lat+u)*r;
    y:=qrand(pos.lat+v, pos.long)*(1.0-r) + qrand(pos.lat, pos.long+u)*r;
  ELSE
    x:=qrand(pos.long, pos.lat);
    y:=qrand(pos.lat, pos.long);    
  END;
  a:=waveamp*one(pixs*(0.5*WAVLEN)/wavelenx);
  wavx:=a*x;
  wavy:=a*y;
END waveang;


PROCEDURE wateralias(wy:REAL; wg-:COL; VAR c:COL);  (* add mirrord image to water colour *)
VAR m:REAL;
BEGIN
  m:=(0.7-wy*0.5);
  IF m<0.1 THEN m:=0.1 END;
  m:=m*waterglance;
  c.r:=c.r*m + wg.r;
  c.g:=c.g*m + wg.g;
  c.b:=c.b*m + wg.b;
END wateralias;

PROCEDURE checkwater(VAR atwa:CARD8; slant:REAL); (* water if flat enough *)
VAR h,ws:REAL;
    w:CARD8;
BEGIN
  ws:=ABS(waterslant);               (* >0 check attribute for flatness else generate water *)
  IF ws<1.0 THEN 
    h:=slant-ws;                                                   (* <0 is not flat enough *)   
    w:=TRUNC(one((slant-ws)/(1.0-ws))*100.0);
    IF waterslant<0.0 THEN atwa:=w ELSIF atwa>w THEN atwa:=w END;  (* make or remove water *)
  END;
END checkwater;


PROCEDURE storepix(image:pIMAGE; x,y:CARDINAL; c:COL; alias:REAL; heaven:BOOLEAN);
VAR h:COL;
BEGIN
--IF heaven THEN INC(debugyi) ELSE INC(debugxi) END;
  IF NOT alticolour THEN c.r:=c.r*LUMMUL; c.g:=c.g*LUMMUL; c.b:=c.b*LUMMUL END;

  IF c.r>30000.0 THEN c.r:=30000.0 END;
  IF c.g>30000.0 THEN c.g:=30000.0 END;
  IF c.b>30000.0 THEN c.b:=30000.0 END;
  WITH image^[x][y] DO
    r:=truncc(c.r); g:=truncc(c.g); b:=truncc(c.b) DIV 2*2+ORD(heaven);
  END;

  IF NOT alticolour & (y>0) THEN
    WITH image^[x][y-1] DO
      IF alias<1.0 THEN
        h.r:=VAL(REAL,r);
        h.g:=VAL(REAL,g);
        h.b:=VAL(REAL,b);
        colblend(alias, h, c, c);    (* antialiasing *)
      END;
      r:=truncc(c.r); g:=truncc(c.g); b:=truncc(c.b) DIV 2*2+ORD(ODD(b) OR heaven);
    END;
  END;
END storepix;

PROCEDURE updatepoirect(pos-, eye:POSITIONL); (* build rectangular around image positions to preselect pois *)
VAR p:POSITION;
BEGIN
  p.lat :=pos.lat -eye.lat;
  p.long:=pos.long-eye.long;
  IF p.lat <poirect0.lat  THEN poirect0.lat :=p.lat END;
  IF p.long<poirect0.long THEN poirect0.long:=p.long END;
  IF p.lat >poirect1.lat  THEN poirect1.lat :=p.lat END;
  IF p.long>poirect1.long THEN poirect1.long:=p.long END;
END updatepoirect;

(*
PROCEDURE wrxyz(x,y,z:LONGREAL);
VAR lat, long, alt:LONGREAL;
BEGIN
  wgs84rl(x, y, z, lat, long, alt);
  WrStr("[");WrFixed(lat/RAD, 6, 1);WrFixed(long/RAD, 6, 12);WrFixed(alt, 1, 8); WrStr("] ");
END wrxyz;
*)

PROCEDURE Panofind(csvstep:CARDINAL; panpar-:PANOWIN; VAR res:REAL; VAR pos:POSITIONL);

CONST FULLUM=1000.0;   (* full bright in image *)
      MAXHP=0.05;      (* vertical spatial luminance highpass *)
      DDIST=50.0;
      LIGHT90=0.35;    (* luminance at 90 deg sun *)
      SELFSHADOW=0.6;  (* dark/light of surface 90/180 deg to sun *)

VAR xi, yi:CARDINAL;

    void:INTEGER;

    maxdist,
    azid, eled, ele0, slong, clong, xn,yn,zn, azi, suneyeaz,
    dd, dm, ddm, wx, wy, slat, clat, refrac, x0,y0,z0, x1,y1,z1, x0m,y0m,z0m :LONGREAL;

    snowdense, pixsize, dlum, lastlum, lum, llum, tree, altf, sloped,
    space, slant, foga, wavelengthx, wavelengthy, waterdist,
    oldh, space1, watershape, xsun, sunarea, dropsun, waterh, pixsizewater, waterdropsun,
    wooddens, snowdens, xsunm, sunlum, wavx, wavy, coastline:REAL;

    pos0, pos1, waterpos:POSITIONL;

    heaven, firsty, mirrord, heavenc, hitwater, dryrun:BOOLEAN;

    ss,ss1:ARRAY[0..100] OF CHAR; 
    atwa, atwo, attu:CARD8;
    sunrotvec:ROTVEC;
    nvec:NORMVEC;
    hc, boc, watergroundcol, snowc, suncol, suncolair, col, aliascol, citylamp:COL;

--    camvec:LVEC;
--    glance:GLANCESET;
    meta:METAINFO;
    pmeta:pMETAINFO;
    istree, odd:BOOLEAN;
    treehints:TREEHINTS;
    ddstack:DDSTACK;
BEGIN
  FILL(ADR(treehints), 0C, SIZE(treehints));

  pmeta:=ADR(meta);

  IF videolat<=90.0 THEN coslat:=cos(videolat*RAD) ELSE coslat:=cos(panpar.eye.lat) END;
  treerasterlat:=treeraster;
  treerasterlong:=treerasterlat*coslat;
  maxdist:=distancel(panpar.eye, panpar.horizon)*1000.0;
  IF maxdist<50.0 THEN RETURN END;                           (* horizon too near *)

  wavelengthx:=wavelength*WAVLEN;
  wavelengthy:=wavelengthx*coslat;
  azieletorotvec(270*RAD - panpar.sunazi, panpar.sunele, sunrotvec);

--azieletovec(-panpar.sunazi, panpar.sunele, glance.sunvec);

  wgs84sl(panpar.eye.lat, panpar.eye.long, panpar.eyealt, x0,y0,z0);

  azi:=azimuthl(panpar.eye, panpar.horizon)*RAD;
--WrFixed(panpar.eye.lat*180.0/PI, 1, 7);WrStrLn(" lat");
  azid:=panpar.angle*RAD/FLOAT(HIGH(panpar.image^)+1);                (* h rad per pixel *)
  IF panpar.flatscreen<>eCURVED THEN azid:=azid*satan(panpar.angle*(RAD*0.5))/(panpar.angle*(RAD*0.5)) END; (* make same hor. sight angle as on courved film *)
--WrFixed(azid*180.0/PI, 2, 7);WrStrLn(" deg/pixel");
  eled:=azid/panpar.yzoom;
  ele0:=panpar.angle*0.5/panpar.yzoom*FLOAT(HIGH(panpar.image^[0])+1)/FLOAT(HIGH(panpar.image^)+1)*RAD;  

  slat:=sin(-panpar.eye.lat);
  clat:=cos(-panpar.eye.lat);
  slong:=sin(-panpar.eye.long);
  clong:=cos(-panpar.eye.long);

  xi:=0;
  refrac:=panpar.refract*REFRACC;
  sunarea:=sunvisability(panpar.sunazi,panpar.sunele,slat,clat,slong,clong, x0,y0,z0, maxdist, refrac);
  IF verb THEN
    WrFixed(azi/RAD, 1, 1); WrStr(", "); WrFixed(panpar.eyealt, 1, 1);
    WrStrLn(" cam azimuth, alt");
    WrFixed(sunarea*100.0, 1, 1); WrStrLn("% of sun visable");
  END;
  teye.lat:=panpar.eye.lat;
  teye.long:=panpar.eye.long;
  showstamdist:=FLOAT(HIGH(panpar.image^))/panpar.angle*60.0;
  suncolour(panpar.sunele, suncol, sunlum);
--WrFixed(suncol.r, 2, 8);WrFixed(suncol.g, 2, 8);WrFixed(suncol.b, 2, 8); WrStrLn("sc");
  REPEAT
    IF dropshadows THEN Shadowmask(xi) END;
    yi:=0;
    dd:=0.0;
    dlum:=0.0;
    lastlum:=0.0;
    heaven:=FALSE;
    firsty:=FALSE;
    mirrord:=FALSE;
    ddstack.startdist:=0.0;
    REPEAT
debugxi:=xi;
debugyi:=yi;
      dryrun:=FALSE;

-----pixel of image to view angles
      IF panpar.flatscreen=eCAMERA THEN
        rotateview(wx, wy, (azid*(LFLOAT(xi) - LFLOAT(HIGH(panpar.image^)+1)*0.5)),
                           (azid*(LFLOAT(yi) - LFLOAT(HIGH(panpar.image^[0])+1)*0.5)),
                           panpar.elevation*RAD);
      ELSE
        wx:=azid*(LFLOAT(xi) - LFLOAT(HIGH(panpar.image^)+1)*0.5);
        wy:=eled*(LFLOAT(yi) - LFLOAT(HIGH(panpar.image^[0])+1)*0.5);
        IF panpar.flatscreen<>eCURVED THEN 
          wx:=atan(wx);
          wy:=atan(wy);
        END;
        wy:=wy + panpar.elevation*RAD;
        IF ABS(wy)>89.0*RAD THEN dryrun:=TRUE END;
      END;

      wx:=wx + azi;
-----pixel of image to view angles

      IF mirrord THEN  (* hit water surface *)
--        wy:=-wy+waveamp*sin(pos.lat*wavelengthx)+waveamp*0.5*sin(pos.lat*wavelengthx*3.7)+waveamp*0.6*sin(pos.long*wavelengthy);
--        wx:= wx+0.5*(waveamp*sin(pos.long*wavelengthx)+0.4*waveamp*sin(pos.long*wavelengthx*2.9)+waveamp*0.6*sin(pos.lat*wavelengthy));
--        wy:=-wy+waveamp*qrand(pos.long*wavelengthy, pos.lat*wavelengthx);
--        wx:= wx+sin(wy)*waveamp*qrand(pos.lat*wavelengthx, pos.long*wavelengthy);
        wx:= wx+sin(wy)*wavy;
        wy:=-wy+wavx;
      END;
--      IF mirrord OR NOT dryrun & ((panpar.flatscreen=eCAMERA) OR (wy>=-89.9999*RAD) & (wy<89.9999*RAD)) THEN

      IF NOT dryrun THEN
        IF NOT heaven THEN

--WrFixed(wx/RAD, 3, 10);WrFixed(wy/RAD, 3, 10); WrStrLn("=wxwy");
          angtoxyz(wx,wy, slat,clat,slong,clong, xn,yn,zn);
          IF mirrord THEN
            dm:=0.0;
--          raytrace(x0m,y0m,z0m, camvec.x,camvec.y,camvec.z, maxdist-dd, -eled, refrac, dm, space1, pos);
            raytrace(x0m,y0m,z0m, xn,yn,zn, maxdist-dd, -eled, refrac, treesize, TRUE, dm, space1, pos, treehints, ddstack, istree, odd, FALSE);
--WrFixed(dm, 1,8); wrxyz(x0m,y0m,z0m);wrxyz(x0m+xn,y0m+yn,z0m+zn);WrFixed(pos.lat/RAD, 6,12);WrFixed(pos.long/RAD, 6,12);WrStrLn("");
            ddm:=dd+dm;
            IF dm>maxdist-dd THEN heaven:=TRUE END;
          ELSE
--          raytrace(x0,y0,z0, camvec.x,camvec.y,camvec.z, maxdist, eled, refrac, dd, space, pos);
            raytrace(x0,y0,z0, xn,yn,zn, maxdist, eled, refrac, treesize, FALSE, dd, space, pos, treehints, ddstack, istree, odd, FALSE);
            ddm:=dd;
            IF dd>maxdist THEN heaven:=TRUE END;
            IF (csvstep<>0) & (xi MOD csvstep=0) (*& (panpar.flatscreen<>eCAMERA)*) THEN xycsv(xi, yi, pos, heaven, firsty) END;
          END;
          IF ddm<1.0 THEN RETURN END;      (* under ground *)
        END;

        lum:=ddm*0.2;
        suneyeaz:=radmod(wx-panpar.sunazi);

(*
        IF NOT mirrord THEN
          llum:=(lum-lastlum)*20.0;
          IF llum>MAXHP*maxdist THEN llum:=MAXHP*maxdist END;
          IF llum>dlum THEN dlum:=llum END;                       (* dlum is highpass peak *) 
          lastlum:=lum;                                           (* lum is distance *)
        END; 
*)
        heavenc:=heaven; 
        IF NOT heaven THEN    
          normvec(pos, pmeta, slant, oldh, nvec, atwa, atwo, attu);
          rotveclight(xsun, nvec, sunrotvec);                     (* cos sunangle on ground *)
          pixsize:=FLOAT(HIGH(panpar.image^))/(ddm*panpar.angle)*(1.0/RAD);    (* pixel/m *)
          IF istree THEN
--            normvectree(treehints[odd]);
--            rotveclight(xsuntree, treehints[odd].nvec, sunrotvec);(* cos sunangle on tree *)
            oldh:=treehints[odd].rooth+treehints[odd].hith;
--            treeignor:=treesize*one(1.0-pixsize*1.5);
          END;
          IF dropshadows & NOT alticolour THEN
            IF xsun>-0.01 THEN                        (* not selfshadow *)
              dropsun:=seesun(pos, oldh+WATERMIRRORGAP, panpar.sunazi, panpar.sunele+SUNSIZE, maxdist, refrac, pixsize);
            ELSE dropsun:=0.0 END;
          ELSE dropsun:=1.0 END;
          IF panpar.sunele<0.0 THEN dropsun:=dropsun+panpar.sunele*15.0 END;  (* night light *)
          shadowcolour(dropsun, suncol, suncolair);
--        shadowcolour(1.0+panpar.sunele*15.0, suncol, suncolair);
          dropsun:=fogsun(dropsun, maxdist, panpar.sunele*RAD, oldh);

          IF alticolour THEN                                      (* only colorized altitude *)
            hc.r:=(xsun+1.0)*(ALTICOLORANG*0.5);
            hc.g:=oldh;
            hc.b:=0.0;
            IF atwa>50 THEN hc.b:=2.0;                            (* store attributes in blue *)
            ELSIF atwo>50 THEN hc.b:=4.0;
            ELSIF attu>50 THEN hc.b:=6.0;
            END;
            IF oldh<alticolmin THEN alticolmin:=oldh END;
            IF oldh>alticolmax THEN alticolmax:=oldh END;

          ELSE
            wooddens:=wooddensity(atwo, oldh, slant);
-- lake ?
            IF NOT mirrord & NOT istree & flatwater THEN checkwater(atwa, slant) ELSE atwa:=0 END;
            IF atwa>0 THEN
              pixsizewater:=pixsize*ABS(sin(wy));
              IF pixsizewater<0.03333 THEN pixsizewater:=0.03333 END;   (* size of srtm is limit for aliasing *)
              waveang(pos, waveamp, wavelengthx, wavelengthy, pixsize, wavx, wavy);  (* add water waves *)
              coastline:=0.5+FLOAT(atwa)-50.0+qrand(pos.long*3000000.0, pos.lat*3000000.0)*2.0;
            ELSE coastline:=-50.0; wavx:=0.0; wavy:=0.0 END;
            hitwater:=coastline>(wavx+wavy)*50.0;

            IF istree THEN
              treecol(pos, treehints[odd], sunrotvec, pixsize, sunlum, wooddens,
              one(1.0-(snowalt-oldh)/slowaltslant), suneyeaz, xsun, wy, atwo, suncolair, hc);
            ELSE
              IF NOT hitwater THEN
                meadowcol(pos, slant, oldh, pixsize, xsun, xsunm, col, coastline);       (* meadow with snow destruction *)
                rockcol(pos, slant, oldh, pixsize, xsun, atwo, istree, (*glance,*) xsunm, col); (* green/rock *)
                citycol(attu, pos, pixsize, wy, col, citylamp);                         (* city *)
                snowdense:=snowdensity(pos, slant, oldh, pixsize, treehints[odd], atwo);
--                snowcol(snowc, suncolair,  xsun-0.05, sunlum);
                snowcol(snowc, suncolair, snownoise(pos, xsun, pixsize), sunlum);
                botaniccol(boc, suncolair, xsunm, sunlum);                     (* apply sun *)
                citylight(hc, col, boc, citylamp, attu);
                colblend(snowdense, hc, snowc, hc);                           (* add snow *)
                darkundertree(atwo, wooddens, snowdense, pixsize, treehints[odd], hc); (* dark under tree *)
              ELSE
                waterpos:=pos;
                waterh:=oldh+WATERMIRRORGAP;
                waterdropsun:=dropsun;
                waterdist:=one(ddm/maxdist)+0.01;
                wgs84sl(pos.lat, pos.long, waterh, x0m,y0m,z0m);
                dm:=0.0;
                mirrord:=waterglance>0.0;
                IF waterslant<0.0 THEN watershape:=1.0 ELSE watershape:=VAL(REAL, atwa)*0.01 END;
                botaniccol(boc, suncolair, sunrotvec.sy, sunlum);                     (* apply sun *)
                underwatercol(pos, pixsize, watershape, snowalt-oldh, wy, wavx, wavy, boc, watergroundcol);
                hc:=watergroundcol;          (* if mirror off *)
              END;
            END;
            IF NOT hitwater THEN
              IF mirrord THEN
                sundisk(suneyeaz, wy, panpar.sunazi, panpar.sunele, one(ddm/maxdist)+0.01,
                        sunarea, maxdist, FALSE, suncol, refrac, hc);
                wateralias(wy, watergroundcol, hc);
                waterglitter(suneyeaz, wy, panpar.sunazi, panpar.sunele, pixsize, maxdist,
                             suncol, waterpos, refrac, waterh, waterdropsun, hc);
                mirrord:=FALSE;
              ELSE waterdist:=one(ddm/maxdist)+0.01 END;
              sundisk(suneyeaz, wy, panpar.sunazi, panpar.sunele, waterdist,
                      sunarea, maxdist, FALSE, suncol, refrac, hc);
            END;

          END;
        ELSE                                                                 (* heaven *)
          IF mirrord THEN 
            sundisk(suneyeaz, wy, panpar.sunazi, panpar.sunele, 0.0,
                    sunarea, maxdist, FALSE, suncol, refrac, hc);
            wateralias(wy, watergroundcol, hc);
            waterglitter(suneyeaz, wy, panpar.sunazi, panpar.sunele, pixsize, maxdist,
                         suncol, waterpos, refrac, waterh, waterdropsun, hc);
            dlum:=0.0;
--WrStrLn("5>");
            heaven:=FALSE;
          ELSE
            waterdist:=0.0;
            sundisk(suneyeaz, wy, panpar.sunazi, panpar.sunele, waterdist, sunarea, maxdist,
                    NOT mirrord, suncol, refrac, hc);
          END; 
          mirrord:=FALSE;
--          heaven:=FALSE;
        END;

IF mirrord & hitwater THEN aliascol:=hc END;

        IF NOT mirrord THEN

(*
IF NOT alticolour & (aliassum>0.0) THEN
--colblend(aliassum, aliascol, hc, hc);
aliassum:=0.0;
--hc.r:=hc.r*aliassum;
--hc.g:=hc.g*aliassum;
--hc.b:=hc.b*aliassum;
END;
*)
          hitwater:=FALSE;
          storepix(panpar.image, xi,yi, hc, space, heavenc);
          IF (panpar.flatscreen=eCAMERA) & (poifiles<>NIL)  THEN updatepoirect(pos, panpar.eye) END;
          dlum:=dlum*(1.0-ddm/maxdist*0.2-120.0/LFLOAT(HIGH(panpar.image^[0])));                                         (* fade out distance highpass peak *)
        END;
      END; 
      IF NOT mirrord THEN
        IF ddstack.startdist<>0.0 THEN                       (* jump back to transparency entry *)
          dd:=ddstack.startdist;
          ddstack.startdist:=0.0;
          heaven:=FALSE;
        END;
        INC(yi);
        firsty:=TRUE;
      END;

    UNTIL yi>HIGH(panpar.image^[0]);




IF verb THEN WrInt(xi, 5); WrStr(CHR(13)); Flush END;
    INC(xi);
  UNTIL xi>HIGH(panpar.image^);
END Panofind;

PROCEDURE Panorama1(panpar-:PANOWIN);
VAR res:REAL;
    pos:POSITIONL;
BEGIN
  Panofind(rastercsv, panpar, res, pos);
END Panorama1;


--- pois


PROCEDURE drawpoiicon(xi, yi:REAL; panpar:PANOWIN; sym-:SYMBOL; xs,ys,cr,cg,cb,calpha:CARDINAL):BOOLEAN;
VAR x,y,iy,ix,iyr,a,ar,lum:CARDINAL;
    ok:BOOLEAN;
BEGIN
  ok:=FALSE;
  IF (xi<=0.0) OR (yi<=0.0) THEN RETURN FALSE END;
  lum:=panpar.maxlum DIV 256;
  cr:=cr*lum;
  cg:=cg*lum;
  cb:=cb*lum;
  x:=VAL(CARDINAL, xi+0.5);
  IF (x>=xs) & (x+xs<=HIGH(panpar.image^)) THEN
    x:=VAL(CARDINAL, xi+0.5-FLOAT(xs)*0.5);
    y:=VAL(CARDINAL, yi+0.5);
    FOR iy:=0 TO ys-1 DO
      iyr:=ys-iy-1;
      IF iy+y<=HIGH(panpar.image^[0]) THEN
        FOR ix:=0 TO xs-1 DO
          WITH panpar.image^[x+ix][y+iy] DO
            a:=sym[iyr][ix].alpha*calpha DIV 256;
            ar:=255-a;
            r:=(a*VAL(CARDINAL, sym[iyr][ix].r8)*cr DIV 256 + ar*VAL(CARDINAL,r)) DIV 255;
            g:=(a*VAL(CARDINAL, sym[iyr][ix].g8)*cg DIV 256 + ar*VAL(CARDINAL,g)) DIV 255;
            b:=(a*VAL(CARDINAL, sym[iyr][ix].b8)*cb DIV 256 + ar*VAL(CARDINAL,b)) DIV 255;
          END;
        END;
        ok:=TRUE;
      END;
    END;
  END;
  RETURN ok
END drawpoiicon;



PROCEDURE drawpoi(xi, yi:REAL; panpar:PANOWIN; size:CARDINAL; cr,cg,cb,calpha:INTEGER):BOOLEAN;
CONST MH=4096;
VAR x,y,iy,ix,lum:CARDINAL;
    c,s,n:INTEGER;
    a:INTEGER;
    ok:BOOLEAN;
BEGIN
  ok:=FALSE;
  IF (xi<=0.0) OR (yi<=0.0) THEN RETURN FALSE END;
  lum:=panpar.maxlum DIV 1024;
  cr:=VAL(CARDINAL,cr)*lum;
  cg:=VAL(CARDINAL,cg)*lum;
  cb:=VAL(CARDINAL,cb)*lum;

  x:=VAL(CARDINAL, xi+0.5);
  y:=VAL(CARDINAL, yi+0.5);
  
  IF (x>=size) & (x+size<=HIGH(panpar.image^)) THEN
    FOR iy:=0 TO size*2-1 DO
      IF iy+y<=HIGH(panpar.image^[0]) THEN
        FOR ix:=0 TO size*2-1 DO
          a:=ABS(VAL(INTEGER,ix)-VAL(INTEGER,size));
          c:=0;
          a:=VAL(INTEGER, iy)-2*a;
          IF a=-1 THEN c:=256*MH DIV 1024
          ELSIF a=0 THEN c:=768*MH DIV 1024
          ELSIF a>0 THEN c:=1024*MH DIV 1024 END;

          IF c>0 THEN
            s:=0;
            IF c<MH THEN
              IF ix<size THEN
                IF iy<>size*2-1 THEN s:=200*MH DIV 1024 END;
              ELSE s:=-400*MH DIV 1024 END; 
            ELSIF iy=size*2-1 THEN c:=500*MH DIV 1024; s:=-300*MH DIV 1024 END;
            c:=c*calpha DIV 256;
            WITH panpar.image^[x+ix-size][y+iy] DO
              n:=(MH-c)*VAL(INTEGER,r) DIV MH + c*cr DIV 256 + s;
              IF n<0 THEN r:=0 ELSE r:=n END;
              n:=(MH-c)*VAL(INTEGER,g) DIV MH + c*cg DIV 256 + s;
              IF n<0 THEN g:=0 ELSE g:=n END;
              n:=(MH-c)*VAL(INTEGER,b) DIV MH + c*cb DIV 256 + s;
              IF n<0 THEN b:=0 ELSE b:=n END;
            END;
          END;
        END;
        ok:=TRUE;
      END;
    END;
  END;
  RETURN ok
END drawpoi;



PROCEDURE drawvispois(panpar-:PANOWIN);  (* draw stored pois if free place *)
VAR p:pPOISTORE;
    poiy, lum:CARDINAL;
    ok:BOOLEAN;
     
BEGIN
  p:=poistore;
  lum:=panpar.maxlum DIV 2;
  WHILE p<>NIL DO
    WITH p^.pp^ DO
      IF p^.textposx>=0 THEN
        poiy:=ys;
        ok:=FALSE;
        IF xs=0 THEN
          ok:=drawpoi(p^.x, p^.y, panpar, poisize, r,g,b,alpha);                    (* internal icon *)
          poiy:=poisize*2;
        ELSIF p^.thrutree THEN
          IF withwoodicon THEN
            ok:=drawpoiicon(p^.x, p^.y, panpar, woodsymbol, xs, ys, r,g,b,alpha);     (* behind wood icon *)
          END;
        ELSE ok:=drawpoiicon(p^.x, p^.y, panpar, symbol, xs, ys, r,g,b,alpha) END;  (* icon from image *)

--- vertical text to icon
        IF ok & (poilabel>0) THEN
          writestr(panpar.image, p^.textposx, truncc(p^.y)+poiy+2, fonttyp, 1, 5,
                   lum, lum, lum, p^.text);
        END;
      END;

    END;
--- vertical text to icon

    p:=p^.next;
  END;
END drawvispois; 

PROCEDURE drawsortpois(panpar-:PANOWIN);  (* draw stored pois if free place *)
VAR p, pn:pPOISTORE;
    poiy:CARDINAL;
    ok, ln, rn:BOOLEAN;
    tx, left, right, i:INTEGER;
    fonty,
    textx, yu, yo:CARDINAL;
    m:SET32;
BEGIN
  p:=poistore;
  poistore:=NIL;
  fonty:=fontsizey(fonttyp);
  WHILE p<>NIL DO
    WITH p^.pp^ DO
      tx:=VAL(INTEGER, p^.x)-VAL(INTEGER, fonty DIV 2);
      IF tx>=0 THEN
        textx:=strsize(fonttyp, p^.text); 
        poiy:=ys;
        IF xs=0 THEN poiy:=poisize*2 END;
        yu:=(truncc(p^.y)+poiy)*32 DIV (HIGH(panpar.image^[0])+1);
        yo:=((truncc(p^.y)+poiy+2+textx)*32+1) DIV (HIGH(panpar.image^[0])+1); 

        IF yu>31 THEN yu:=31 END;
        IF yo>31 THEN yo:=31 END;
        m:=SET32{};
        FOR i:=yu TO VAL(INTEGER,yo) DO INCL(m, i) END;
--      m:=SET32{yo..yu};
        IF tx>VAL(INTEGER, HIGH(poimap^)) THEN tx:=HIGH(poimap^) END;

        left:=tx-VAL(INTEGER,fonty DIV 3);          (* less activity in lowest pixel rows *)
        IF left<0 THEN left:=0 END;
        right:=tx+VAL(INTEGER,fonty+fonty DIV 2);
        IF right>VAL(INTEGER, HIGH(poimap^)) THEN right:=HIGH(poimap^) END;
        ok:=FALSE;
        rn:=FALSE;
        ln:=FALSE;
        i:=left; 
        LOOP
          IF poimap^[i]*m<>SET32{} THEN                       (* not full free image column *)
            IF ok THEN right:=i-1; rn:=TRUE; EXIT END;

            left:=i+1;
            ln:=TRUE;
          ELSE ok:=TRUE END;
          INC(i);
          IF i>=right THEN EXIT END;
        END;
        IF ln THEN tx:=left;  (* dock left or right *)
        ELSIF rn THEN tx:=imax(tx-VAL(INTEGER,fonty DIV 2), right-VAL(INTEGER,fonty)) END;
        IF (tx>=0) & ((poilabel<=1) OR (right-left>=VAL(INTEGER,fonty))) THEN (* text fits in empty place *) 
          right:=tx+VAL(INTEGER,fonty)-1;
          IF right>VAL(INTEGER, HIGH(poimap^)) THEN right:=HIGH(poimap^) END;
          FOR i:=tx TO right DO poimap^[i]:=poimap^[i]+m END;  (* mark where is text on image *)
        ELSE tx:=-1 END;
        p^.textposx:=tx;
      END;
    END;
    pn:=p^.next;
    p^.next:=poistore;         (* revert direction to overdraw farer icons with neerer *)
    poistore:=p;
    p:=pn;
  END;
END drawsortpois; 

PROCEDURE storepoi(x, y, z:REAL; thrutree:BOOLEAN; text-:ARRAY OF CHAR; pp:pPOI);  (* sort in pois on priority and distance *)
VAR p, po, pn:pPOISTORE;
BEGIN
  p:=poistore;
  po:=NIL;
  WHILE (p<>NIL) & (pp^.prio>=p^.prio) & ((pp^.prio>p^.prio) OR (z>p^.dist)) DO po:=p; p:=p^.next END; 
  ALLOCATE(pn, SIZE(pn^));
  IF pn=NIL THEN Werr("png write out of memory"+LF) END;

  pn^.pp:=pp;
  pn^.x:=x;
  pn^.y:=y;
  pn^.dist:=z;
  pn^.prio:=pp^.prio;
  pn^.thrutree:=thrutree;
  Assign(pn^.text, text);
  pn^.next:=p;
  IF po<>NIL THEN po^.next:=pn ELSE poistore:=pn END;
END storepoi;

PROCEDURE backtrackpixel(pos:POSITIONL; alt:REAL; panpar-:PANOWIN; VAR px, py:REAL;
                    VAR treehints:TREEHINTS; altmode:CARDINAL; VAR rd:LONGREAL);

VAR nne, nnp, aeye, d2, xi, yi, refrac:REAL;
    resoltx, azid, azi, eled, subpix:REAL;

    x0,y0,z0, x1,y1,z1, xn,yn,zn, xnn, ynn, znn:LONGREAL;
    atwa:CARD8;
    rpos:POSITIONL;
    istree, odd, ok:BOOLEAN;
    poiy:CARDINAL;
    xt:INTEGER; 
    ddstack:DDSTACK;
BEGIN
  px:=-1.0;
  py:=-1.0;
  rd:=-1.0;
  IF NOT posvalidl(pos) THEN RETURN END;
  refrac:=panpar.refract*REFRACC;
  nne:=getsrtmlong(panpar.eye.lat, panpar.eye.long, 0, TRUE, resoltx, atwa, NIL);
  IF defaultocean & (nne>NOALT) THEN nne:=0.0 END;
  nnp:=getsrtmlong(pos.lat, pos.long, 0, TRUE, resoltx, atwa, NIL);
  IF defaultocean & (nnp>NOALT) THEN nnp:=0.0 END;
  IF altmode=0 THEN alt:=nnp;                        (* ignore track alt *) 
  ELSIF altmode=1 THEN alt:=max(alt, nnp);           (* minimum alt is from srtm *)  
  ELSIF altmode=2 THEN
    IF alt<nnp THEN RETURN END;                      (* dont fly under ground *)

  ELSIF altmode=3  (* pp^.altIsOG*) THEN alt:=nnp + max(alt, minpoialt);
  ELSE
    nnp:=nnp + minpoialt;
    alt:=max(nnp,alt);                               (* set file defined altitude *)
  END;

  aeye:=panpar.eyealt;
  wgs84sl(panpar.eye.lat, panpar.eye.long, aeye, x0,y0,z0);
  wgs84sl(pos.lat, pos.long, alt, x1,y1,z1);
  d2:=sqrt(sqr(x1-x0) + sqr(y1-y0) + sqr(z1-z0));
  wgs84sl(pos.lat, pos.long, (alt+refrac*d2*d2), x1,y1,z1);           (* apply refraction *)
  d2:=sqrt(sqr(x1-x0) + sqr(y1-y0) + sqr(z1-z0));

  xn:=(x1-x0)/d2;                                                     (* normvector *)
  yn:=(y1-y0)/d2;
  zn:=(z1-z0)/d2;
  xnn:=xn;
  ynn:=yn;
  znn:=zn;
  xn:=-xn;
  rotvector(yn,xn, cos(panpar.eye.long),sin(panpar.eye.long));
  rotvector(zn,xn, cos(panpar.eye.lat), sin(panpar.eye.lat));
  azi:=azimuthl(panpar.eye, panpar.horizon)*RAD;
  rotvector(zn,yn, cos(azi), sin(azi));
  IF panpar.flatscreen=eCAMERA THEN
    rotvector(xn,zn, cos(panpar.elevation*RAD), sin(panpar.elevation*RAD));
  END;
  azid:=panpar.angle*RAD/FLOAT(HIGH(panpar.image^)+1);                (* h rad per pixel *)
  IF panpar.flatscreen<>eCURVED THEN azid:=azid*satan(panpar.angle*(RAD*0.5))/(panpar.angle*(RAD*0.5)) END;
  eled:=azid/panpar.yzoom;                                            (* v rad per pixel *)

  xi:=atang2(zn,yn) (* - azi*);
  IF xi<-PI THEN xi:=xi+2*PI ELSIF xi>PI THEN xi:=xi-2*PI END;
  IF panpar.flatscreen<>eCURVED THEN xi:=satan(xi) END;
  xi:=xi/azid + FLOAT(HIGH(panpar.image^)+1)*0.5;
  IF panpar.flatscreen<>eCAMERA THEN
    yi:=atang2(sqrt(zn*zn + yn*yn),-xn)-panpar.elevation*RAD;
  ELSE
    yi:=atang2(zn,-xn);
  END;
  IF panpar.flatscreen<>eCURVED THEN yi:=satan(yi) END;
  yi:=yi/eled + FLOAT(HIGH(panpar.image^[0])+1)*0.5;
  IF (xi>=0.0) & (xi<=FLOAT(HIGH(panpar.image^)))
   & (yi>=0.0) & (yi<=FLOAT(HIGH(panpar.image^[0]))) THEN
    rd:=0.0;
    FILL(ADR(treehints), 0C, SIZE(treehints));
    odd:=FALSE;
    treehints[FALSE].treepass:=TRUE;
    treehints[FALSE].hith:=treesize;
    ddstack.startdist:=0.0;
    raytrace(x0,y0,z0, xnn,ynn,znn, d2, 0.0, refrac, treesize, FALSE, rd, subpix, rpos, treehints,
           ddstack, istree, odd, FALSE);
--  IF (rd+30.0<d2) OR (altmode=2) & (alt<nnp) THEN rd:=-1.0 END;                                       (* point hidden *)
    IF rd+30.0<d2 THEN rd:=-1.0 END;                                       (* point hidden *)
  END;  
  px:=xi;
  py:=yi;
END backtrackpixel;


PROCEDURE posinimage(panpar-:PANOWIN; pos:POSITIONL):BOOLEAN;
VAR d:LONGREAL;
    az:REAL;
BEGIN
  IF panpar.flatscreen=eCAMERA THEN    (* in rectangular *)
    d:=pos.lat-panpar.eye.lat;
    IF (d<poirect0.lat) OR (d>poirect1.lat) THEN RETURN FALSE END;

    d:=pos.long-panpar.eye.long;
    IF (d<poirect0.long) OR (d>poirect1.long) THEN RETURN FALSE END;

  ELSE                                 (* in circle segment *)
    d:=distancel(panpar.eye, pos);
    IF (d>panpar.dist) OR (d<0.05) THEN RETURN FALSE END;

    az:=(azimuthl(panpar.eye, panpar.horizon)-azimuthl(panpar.eye, pos))*RAD;
    IF az<-PI THEN az:=az+2.0*PI ELSIF az>PI THEN az:=az-2.0*PI END;
    IF ABS(az)>panpar.hseg THEN RETURN FALSE END;

  END;
  RETURN TRUE
END posinimage;

PROCEDURE setpoi(panpar-:PANOWIN; posdeg-:POSITION; alti:REAL; text-, ext-:ARRAY OF CHAR; pp:pPOI);
VAR alt, resol, og:REAL;
    pos:POSITIONL;
    atwa:CARD8;
    px, py:REAL;
    treehints:TREEHINTS;
    dist:LONGREAL;
BEGIN
  pos.lat :=posdeg.lat*RAD;
  pos.long:=posdeg.long*RAD;

--WrFixed(posdeg.lat,2,9); WrStr(" ");
--WrFixed(posdeg.long,2,9); WrStr(" ");

  IF NOT posinimage(panpar, pos) THEN RETURN END;

  backtrackpixel(pos, alti, panpar, px, py, treehints, 4-ORD(pp^.altIsOG), dist);
  IF dist>=0.0 THEN storepoi(px, py, dist, treehints[FALSE].hith<treesize, text, pp) END;

  IF verb THEN

WrFixed(distancel(panpar.eye, pos),3,9); WrStr(" ");
WrFixed(azimuthl(panpar.eye, pos),2,9); WrStr(" ");
IF alti>-1000.0 THEN WrFixed(alti, 0,9); WrStr(" ") END;
WrStr(text);
IF dist<0.0 THEN
  WrStr(" #");
ELSIF treehints[FALSE].hith<treesize THEN
  WrStr(" (");
  WrFixed(treehints[FALSE].hith, 1,1);
  WrStr("m/");
  WrFixed(treehints[FALSE].passdist*0.001, 3,1);
  WrStr("km)");
END;
WrStrLn("");

  END;

  IF dist>=0.0 THEN 
    alt:=getsrtmlong(pos.lat, pos.long, 0, TRUE, resol, atwa, NIL);
    IF defaultocean & (alt>NOALT) THEN alt:=0.0 END;
    IF (alti>=10000.0) OR (alti<alt) THEN alti:=alt END;
    wrcsv(truncc(px+0.5), truncc(py+0.5), posdeg, alti, text, pp^.iconhint, ext);
  END;

END setpoi;



PROCEDURE rdmountain(panpar-:PANOWIN; fn:ARRAY OF CHAR; pp:pPOI); (* import csv file with mountain name, pos, altitude *)
CONST POIINFOSIZE=4096;
VAR p,len,r, sl, si:INTEGER;
    fd:File;
    pos:POSITION;
    alt:REAL;
    b:ARRAY[0..4095] OF CHAR;
    s:ARRAY[0..1023] OF CHAR;
    text:ARRAY[0..POIINFOSIZE] OF CHAR;
    com, lat, long, name, alti, ext:ARRAY[0..256] OF CHAR;


  PROCEDURE getch():CHAR;
  BEGIN
    IF p>=len THEN
      len:=RdBin(fd, b, SIZE(b));
      IF len<=0 THEN RETURN 0C END;

      p:=0; 
    END;
    INC(p);
    RETURN b[p-1]
  END getch;

  PROCEDURE getword(VAR s:ARRAY OF CHAR):INTEGER;
  VAR i:CARDINAL;
      inqu:BOOLEAN;
  BEGIN
    i:=0;
    inqu:=FALSE;
    LOOP
      s[i]:=getch();
      IF s[i]=0C THEN RETURN -1 END;

      IF s[i]=LF THEN s[i]:=0C; RETURN 0 END;

      IF NOT inqu & (s[i]=",") THEN s[i]:=0C; RETURN 1 END;

      IF s[i]='"' THEN inqu:=NOT inqu
      ELSIF (i<HIGH(s)) & (s[i]>=" ") & (s[i]<200C)THEN INC(i) END;
    END;
  END getword; 


BEGIN
  fd:=OpenRead(fn);
  IF NOT FdValid(fd) THEN Error("poi-file not found") END;

  p:=0;
  len:=0;
  LOOP          
    com[0]:=0C;
    name[0]:=0C;
    text[0]:=0C;
    lat[0]:=0C;
    long[0]:=0C;
    alti[0]:=0C;
    ext[0]:=0C;
    r:=getword(com);
    IF r>0 THEN
      r:=getword(name);
      IF r>0 THEN
        r:=getword(text);
        IF r>0 THEN
          r:=getword(lat);
          IF r>0 THEN
            r:=getword(long);
            IF r>0 THEN
              r:=getword(alti);
              IF r>0 THEN
                r:=getword(ext);
                WHILE r>0 DO r:=getword(s) END;
              END;
            END;
          END;
        END;
      END;
    END;
    IF r<0 THEN EXIT END;

    IF (alti[0]=0C) OR NOT StrToFix(alt, alti) THEN alt:=-30000.0 END;

    IF (com[0]<>"#") & (name[0]<>0C) & StrToFix(pos.lat, lat)
    & StrToFix(pos.long, long) & posvalid(pos) THEN
      setpoi(panpar, pos, alt, name, ext, pp);
    END;
  END;
  Close(fd);
END rdmountain;

PROCEDURE rdmountains(panpar-:PANOWIN);
VAR pp:pPOI;
    i:CARDINAL;
BEGIN
  poistore:=NIL;
  NEW(poimap, HIGH(panpar.image^)+1);
  IF poimap=NIL THEN Error("poimap, out of memory") END;

  FOR i:=0 TO HIGH(poimap^) DO poimap^[i]:=SET32{} END;
  pp:=poifiles;
  WHILE pp<>NIL DO
    pp^.xs:=0;
    pp^.ys:=0;
    IF pp^.iconfn[0]<>0C     THEN readsymbol(pp^.symbol, pp^.xs, pp^.ys, pp^.iconfn) END;
    IF pp^.iconwoodfn[0]<>0C THEN
      readsymbol(pp^.woodsymbol, pp^.xs, pp^.ys, pp^.iconwoodfn);
      pp^.withwoodicon:=TRUE;
    ELSE pp^.withwoodicon:=FALSE END;
    rdmountain(panpar, pp^.poifn, pp);

    pp:=pp^.next;
  END;
  drawsortpois(panpar);
  drawvispois(panpar);
END rdmountains;


---- pois

---- tracks




PROCEDURE interpolway(w:REAL; wp:pWAYPOINTS; cnt:CARDINAL; VAR pos:POSITIONL; VAR alt:LONGREAL):BOOLEAN; (* 0..1 for whole track *)
VAR wi:CARDINAL;
    v, x0,y0,z0,x1,y1,z1, a:LONGREAL;   
    ok, noalt:BOOLEAN;
BEGIN
  IF (w<0.0) OR (w>1.0) OR (cnt<2) THEN RETURN FALSE END;

  ok:=FALSE;
  w:=w*FLOAT(cnt-1);
  wi:=truncc(w);
  v:=w-FLOAT(wi);
  WHILE (wp<>NIL) & (wi>HIGH(wp^.wps)) DO DEC(wi, HIGH(wp^.wps)+1); wp:=wp^.next END;
  IF (wp<>NIL) & (wi<cnt) THEN
(*
    x0:=wp^.wps[wi].pos.lat;
    y0:=wp^.wps[wi].pos.long;
    a0:=wp^.wps[wi].alt;
*)
    a:=wp^.wps[wi].alt;
    noalt:=a<-1000.0;
    IF noalt THEN a:=0.0 END;
    wgs84sl(wp^.wps[wi].pos.lat, wp^.wps[wi].pos.long, a, x0,y0,z0);  (* interpolate in wgs84 space for long radio link segments *)
    INC(wi);
    WHILE (wp<>NIL) & (wi>HIGH(wp^.wps)) DO DEC(wi, HIGH(wp^.wps)+1); wp:=wp^.next END;
    IF (wp<>NIL) & (wi<cnt) THEN
(*
      x1:=wp^.wps[wi].pos.lat;
      y1:=wp^.wps[wi].pos.long;
      a1:=wp^.wps[wi].alt;
      pos.lat :=(x0*(1.0-v) + x1*v)*RAD;
      pos.long:=(y0*(1.0-v) + y1*v)*RAD;
      alt :=a0*(1.0-v) + a1*v;
*)
      a:=wp^.wps[wi].alt;
      noalt:=noalt OR (a<-1000.0);
      IF noalt THEN a:=0.0 END;
      wgs84sl(wp^.wps[wi].pos.lat, wp^.wps[wi].pos.long, a, x1,y1,z1);
      wgs84rl(x0*(1.0-v) + x1*v, y0*(1.0-v) + y1*v, z0*(1.0-v) + z1*v, pos.lat, pos.long, alt);
      IF noalt THEN alt:=-30000.0 END;
      ok:=TRUE;
    END;
  END;
  RETURN ok
END interpolway;

PROCEDURE wppixel(x,y:CARDINAL; rr,gg,bb:CARDINAL; alpha:REAL; panpar-:PANOWIN);
VAR a:CARDINAL;
    c, d:CARDINAL;
BEGIN
  WITH panpar.image^[x][y] DO
    a:=truncc(alpha*(64.0*256.0));
    IF a>64*256 THEN a:=64*256 END;
    d:=truncc(256.0-alpha*64.0);            (* darken original *)
    c:=ASH(VAL(CARDINAL,r)*d + rr*a,-8); IF c>32000 THEN c:=32000 END; r:=c;
    c:=ASH(VAL(CARDINAL,g)*d + gg*a,-8); IF c>32000 THEN c:=32000 END; g:=c;
    c:=ASH(VAL(CARDINAL,b)*d + bb*a,-8); IF c>32000 THEN c:=32000 END; b:=c;
  END;
END wppixel;

PROCEDURE bigwp(px,py:REAL; rr,gg,bb:CARDINAL; r:REAL; panpar-:PANOWIN);
VAR x,y, xi, yi,q:INTEGER;
    b,w:REAL;
BEGIN
  x:=truncc(px);
  y:=truncc(py);
  w:=3.0/(r+2.0); 
  q:=truncc(r+2.0);
  FOR yi:=y-q TO y+q DO
    FOR xi:=x-q TO x+q DO
      b:=1.0-sqrt(sqr(FLOAT(xi)-px) + sqr(FLOAT(yi)-py))/(0.5+r);
      IF b>1.0 THEN b:=1.0 END;
      IF (b>0.0) & (xi>=0) & (yi>=0) & (xi<=VAL(INTEGER,HIGH(panpar.image^))) & (yi<=VAL(INTEGER,HIGH(panpar.image^[0]))) THEN
        wppixel(xi,yi, rr,gg,bb, b*w, panpar);
      END;
    END;
  END; 
END bigwp;


PROCEDURE drawtrack(waypoints:pWAYPOINTS; tf:pTRACK; panpar-:PANOWIN);
CONST FRAME=20.0;
VAR wi, wo, d, dp, oldd, oldx, oldy:REAL;
    alt, vis:LONGREAL;
    pos:POSITIONL; 
    px, py, xo, yo, os, limw, pixsize:REAL;
    xf, yf:REAL;
    treehints:TREEHINTS;
    x,y, nok, mok, altm, last:CARDINAL;
    inf, lasrt:BOOLEAN;
    lastdir:INTEGER;
BEGIN
  d:=0.00001;
  wi:=0.0;
  nok:=0;
  mok:=5;
  inf:=FALSE;
  last:=0;
  oldd:=0.0;
  LOOP
    IF wi<0.999999 THEN last:=0 END;
    IF (last<10) & (wi>=1.0) THEN wi:=0.9999999; INC(last) END;
    IF NOT interpolway(wi, waypoints, tf^.cnt, pos, alt) THEN EXIT END;

    altm:=tf^.altmode;
    IF altm=3 THEN altm:=2 END;
    backtrackpixel(pos, alt, panpar, px, py, treehints, altm, vis);
    pixsize:=FLOAT(tf^.width)*min(2.0,0.5+panpar.dist/max(vis,1.0)*50.0);
    limw:=min(3.0, max(1.0, pixsize-1.0));
    os:=max(-px, 0.0);
    os:=max(os, -py);
    os:=max(os, px-FLOAT(HIGH(panpar.image^)-1));
    os:=max(os, py-FLOAT(HIGH(panpar.image^[0])-1));
    os:=min(os, 1000.0);                                (* how many pixels outside image *)
    IF (os<1.0) & NOT inf THEN d:=0.00001 END;          (* immediate fine steps from out to inside image *)

    IF (wi<>0.0) & (nok<=mok) THEN                      (* dont retry forever *)
      IF inf OR (os<1.0) THEN 
        dp:=sqrt(sqr(px-xo)+sqr(py-yo));
        IF dp<0.8*limw THEN                             (* hold waypoint steps in tolerance *)
          INC(nok);
          oldd:=d;
          oldx:=px;
          oldy:=py; 
          d:=d*1.2;
        ELSIF dp>1.3*limw THEN
          INC(nok);
          d:=d*0.9;
          IF (dp>5.0*limw) (* OR (oldd<>0.0)*) THEN
            d:=0.000001;
            wo:=wo+oldd;
            xo:=oldx;
            yo:=oldy;
          END;
          oldd:=0.0;
        ELSE
          nok:=0;
          oldd:=0.0;
        END;
        IF nok<>0 THEN                                  (* retry different step from old start *)
          wi:=wo;
          px:=xo;
          py:=yo;
        END;
      ELSE d:=0.0001+os*0.000001 END;                   (* larger steps far outside image *)

    ELSIF nok>mok THEN nok:=0 END;
    IF nok=0 THEN
      xo:=px;                                           (* store last good step start *)
      yo:=py;
      wo:=wi;
      mok:=5;
      IF (vis>=0.0) & (px>=0.0) & (py>=0.0) & (px<FLOAT(HIGH(panpar.image^))) & (py<FLOAT(HIGH(panpar.image^[0]))) THEN  (* destination not hidden *)
        mok:=50;
(*
        x:=truncc(px);
        xf:=px-FLOAT(x);
        y:=truncc(py);
        yf:=py-FLOAT(y);
        wppixel(x+1, y+1, tf^.r,tf^.g,tf^.b, xf*yf, panpar);            (* draw antialised pixel *)
        wppixel(x  , y+1, tf^.r,tf^.g,tf^.b, (1.0-xf)*yf, panpar);
        wppixel(x+1, y,   tf^.r,tf^.g,tf^.b, xf*(1.0-yf), panpar);
        wppixel(x  , y,   tf^.r,tf^.g,tf^.b, (1.0-xf)*(1.0-yf), panpar);
*)

        bigwp(px,py, tf^.r,tf^.g,tf^.b, pixsize, panpar);
      END;
    END;
    inf:=os<1.0;
    wi:=wi+d;
  END;
END drawtrack;


PROCEDURE rdtracks(panpar-:PANOWIN);
VAR tf:pTRACK;
    b:ARRAY[0..4095] OF CHAR;
    s:ARRAY[0..1023] OF CHAR;
    lat, long, alti:ARRAY[0..256] OF CHAR;
    waypoints, pt, pto:pWAYPOINTS;
    p, len, wp, subtr, cnt:CARDINAL;
    r,rr:INTEGER;
    bl:BOOLEAN;
    nnp, newalt, resoltx, dst:REAL;
    atwa:CARD8;
    newpos:POSITIONL;
 

  PROCEDURE getch():CHAR;
  BEGIN
    IF p>=len THEN
      len:=RdBin(tf^.fd, b, SIZE(b));
      IF len<=0 THEN RETURN 0C END;

      p:=0;
    END;
    INC(p);
    RETURN b[p-1]
  END getch;

  PROCEDURE getword(VAR s:ARRAY OF CHAR):INTEGER;
  VAR i:CARDINAL;
      inqu:BOOLEAN;
  BEGIN
    i:=0;
    inqu:=FALSE;
    LOOP
      s[i]:=getch();
      IF s[i]=0C THEN RETURN -1 END;

      IF s[i]=LF THEN s[i]:=0C; RETURN 0 END;

      IF NOT inqu & (s[i]=",") THEN s[i]:=0C; RETURN 1 END;

      IF s[i]='"' THEN inqu:=NOT inqu
      ELSIF (i<HIGH(s)) & (s[i]>=" ") & (s[i]<200C) THEN INC(i) END;
    END;
  END getword;


BEGIN
  tf:=trackfiles;
  p:=0;
  len:=0;
  subtr:=0;
  WHILE tf<>NIL DO
    waypoints:=NIL;
    wp:=HIGH(pt^.wps)+1;
    tf^.cnt:=0;
    pto:=NIL;
    bl:=FALSE;
    cnt:=0;
    LOOP
      lat[0]:=0C;
      long[0]:=0C;
      alti[0]:=0C;
      r:=getword(lat);
      IF r>0 THEN
        r:=getword(long);
        IF r>0 THEN
          rr:=getword(alti);
          WHILE rr>0 DO rr:=getword(s) END;
        END;
      ELSIF r=0 THEN bl:=TRUE; EXIT END;         (* blank line in csv is start new track *)
      IF r<0 THEN EXIT END;
      INC(cnt);
      IF lat[0]<>"#" THEN
        IF StrToFixL(newpos.lat, lat)
        & StrToFixL(newpos.long, long) & posvalidl(newpos)
        & (ABS(newpos.lat)<90.0) & (ABS(newpos.long)<=180.0) THEN

          newpos.lat:=newpos.lat*RAD;
          newpos.long:=newpos.long*RAD;

          dst:=distancel(panpar.eye, newpos);
          IF dst<panpar.dist THEN

            IF (alti[0]=0C) OR NOT StrToFix(newalt, alti)
            OR (newalt<-1000.0) THEN newalt:=-30000.0 END;

            IF (newalt>=-1000.0) & (tf^.altmode=3) THEN    (* make over ground *)
              nnp:=getsrtmlong(newpos.lat, newpos.long, 0, TRUE, resoltx, atwa, NIL);
              newalt:=nnp + max(newalt, minpoialt);
            END;

            IF wp>HIGH(pt^.wps) THEN
              NEW(pt);
              IF pt=NIL THEN Error("tracks, out of memory") END;

              wp:=0;
              IF pto=NIL THEN waypoints:=pt ELSE pto^.next:=pt END;
              pto:=pt;
              pt^.next:=NIL; 
            END;

            pt^.wps[wp].pos:=newpos;
            pt^.wps[wp].alt:=newalt;

            INC(wp);
            INC(tf^.cnt);
          ELSIF verb THEN WrStr("#") END;
IF verb THEN
  WrStr(tf^.fn); 
  IF subtr>0 THEN WrStr("("); WrInt(subtr,1); WrStr(")") END; 
  WrFixed(newpos.lat*(1.0/RAD), 5,10);
  WrFixed(newpos.long*(1.0/RAD), 5,11);
  IF newalt>-1000.0 THEN WrFixed(newalt, 1,8);WrStr("m"); END;
  WrStr(" "); WrFixed(dst, 1,1);
  WrStrLn("km");
END;

        ELSIF verb THEN WrStr("---csv error line:"); WrInt(cnt,1); WrStrLn("") END;
      END;
    END;
    drawtrack(waypoints, tf, panpar);
    WHILE waypoints<>NIL DO pt:=waypoints; waypoints:=waypoints^.next; DISPOSE(pt) END; 
    IF NOT bl THEN                    (* else go ahead in same file *)
      tf:=tf^.next;
      p:=0;
      len:=0;
    ELSE INC(subtr) END;
  END;
END rdtracks;


---- tracks

PROCEDURE WrLatLong(panpar-:PANOWIN);
CONST DEGSYM=177C;

VAR y,lum:CARDINAL;
    s,h:ARRAY[0..200] OF CHAR;
BEGIN
  FixToStr(panpar.eye.lat/RAD  , 6, s); Append(s, "/");
  FixToStr(panpar.eye.long/RAD , 6, h); Append(s, h); Append(s, " A=");
  FixToStr(panpar.eyealt       , 2, h); Append(s, h); Append(s, " OG=");
  FixToStr(panpar.eyeog        , 2, h); Append(s, h); Append(s, " d=");
  FixToStr(panpar.dist         , 2, h); Append(s, h); Append(s, " e=");
  FixToStr(panpar.elevation    , 2, h); Append(s, h); Append(s, DEGSYM);
  IF panpar.flatscreen=eCAMERA THEN
    Append(s, " center="); FixToStr(azimuthl(panpar.eye, panpar.horizon), 2, h); Append(s, h); Append(s, DEGSYM);
    Append(s, " angle="); FixToStr(panpar.angle, 2, h); Append(s, h); Append(s, DEGSYM);
  END;
  Append(s, " sun=");
  FixToStr(panpar.sunazi/RAD   , 2, h); Append(s, h); Append(s, "/");
  FixToStr(panpar.sunele/RAD   , 2, h); Append(s, h); Append(s, DEGSYM);
  IF alticolour THEN
    Append(s, " altcol=");
    FixToStr(alticolmin, 0, h); Append(s, h); Append(s, ":");
    FixToStr(alticolmax, 0, h); Append(s, h);
    IF altcolstep>0 THEN Append(s, "/"); IntToStr(altcolstep, 0, h); Append(s, h) END;
    Append(s, "m");
  ELSE
    Append(s, " snow=");
    FixToStr(snowalt, 0, h); Append(s, h); Append(s, "/");
    FixToStr(slowaltslant, 0, h); Append(s, h); Append(s, "m W=");
    FixToStr(waveamp*100.0, 2, h); Append(s, h); Append(s, "/");
    FixToStr(wavelength, 2, h); Append(s, h); Append(s, "/");
    FixToStr(waterslant, 3, h); Append(s, h); Append(s, " T=");
    FixToStr(treesize, 0, h); Append(s, h); Append(s, "/");
    FixToStr(botanicslant, 0, h); Append(s, h); Append(s, "/");
    FixToStr(treerastersize, 1, h); Append(s, h); Append(s, "/");
    FixToStr(zerotreealt, 0, h); Append(s, h); Append(s, "/");
    FixToStr(fulltreealt, 0, h); Append(s, h); Append(s, "m");
  END;
  IF dateh[0]<>0C THEN
    Append(s, " utc:"); Append(s, dateh);
    Append(s, " "); Append(s, timeh);
  END;
  lum:=panpar.maxlum DIV 6;
  y:=1;
  IF panpar.flatscreen<>eCAMERA THEN y:=fontsizey(fonttypfoot) END;
  writestr(panpar.image,5, y, fonttypfoot, 0, 10, lum,lum,lum, s);
END WrLatLong;


PROCEDURE drawscale(panpar-:PANOWIN; px, deg:CARDINAL);
VAR y, lum:CARDINAL;
    s:ARRAY[0..100] OF CHAR;
BEGIN
  lum:=panpar.maxlum DIV 6;

  IF (px>2) & (px+2<=HIGH(panpar.image^)) THEN
    FOR y:=0 TO scalesize DIV 10 DO
      WITH panpar.image^[px][y] DO
        r:=r+lum;
        g:=g+lum;
      END; 
    END;
    IF px<=HIGH(panpar.image^) THEN
      IF deg=0 THEN s:="N"
      ELSIF deg=270 THEN s:="W" ELSE IntToStr(deg, 0, s) END; 
      writestr(panpar.image,px+2, 0, fonttypfoot, 0, 5, lum,lum,lum, s);
    END;
  END;
END drawscale;


PROCEDURE scale(panpar-:PANOWIN);
VAR a,b,d,w,i:CARDINAL;
    azi, azid, wx:REAL; 
BEGIN
  azid:=panpar.angle/FLOAT(HIGH(panpar.image^)+1);              (* deg / pixel *)
  d:=TRUNC(FLOAT(scalesize)*azid);
  a:=1;
  b:=1;
  WHILE a*b<d DO
    IF a=1 THEN a:=2 ELSIF a=2 THEN a:=5 ELSIF a=5 THEN a:=1; b:=b*10 END;
  END;
  a:=a*b;                                                       (* 1 2 5 10 deg / step *) 
  IF a>90 THEN a:=180 ELSIF a>60 THEN a:=90 ELSIF a>20 THEN a:=60 END;
--WrInt(a, 10); WrStrLn(" deg/ step");
  azi:=azimuthl(panpar.eye, panpar.horizon);
  wx:=azi-panpar.angle*0.5;
  IF wx<0.0 THEN wx:=wx+360.0 END;
  w:=VAL(CARDINAL, wx) DIV a*a;
  FOR i:=0 TO (VAL(CARDINAL,panpar.angle)+a) DIV a DO
--WrInt(w, 4);
    wx:=(FLOAT(w)-azi)*RAD;
    IF wx<-PI THEN wx:=wx+2*PI ELSIF wx>PI THEN wx:=wx-2*PI END;
    IF panpar.flatscreen<>eCURVED THEN wx:=satan(wx) END;
    wx:=(wx*(1.0/RAD)+panpar.angle*0.5)/azid;
    IF wx>0.0 THEN drawscale(panpar, VAL(CARDINAL, wx+0.5), w) END;
--WrFixed(wx,1, 7);

    INC(w, a);
    IF w>=360 THEN DEC(w, 360) END;
  END;
  WrLatLong(panpar);
--WrStrLn("");
END scale;

PROCEDURE colourize(img:pIMAGE; amin, amax:REAL; atab-:ALTCOLTAB);  (* make altitude to colour *)
CONST LUM=6;
      SHADOWLUM=ALTICOLORANG DIV 60000;
      TAGS=6;
      MAXLUT=1024;

 PROCEDURE intpol(ca,cb:CARDINAL; f:REAL):REAL;
 BEGIN RETURN (LUM/255.0)*(FLOAT(ca)*(1.0-f)+FLOAT(cb)*f) END intpol; 

VAR amul, a, f, K:REAL;
    x, y, c:CARDINAL;
    h, hs:INTEGER;
    lut:ARRAY[0..MAXLUT-1] OF RECORD rr,gg,bb:REAL END;
    tlut:ARRAY[0..TAGS-1] OF RECORD on:BOOLEAN; rr,gg,bb:REAL END;

BEGIN
  IF atab.len=0 THEN                                (* no file make default lut *)
    K:=LUM/FLOAT(MAXLUT);
    FOR c:=0 TO HIGH(lut) DO
      WITH lut[c] DO
        IF c<400 THEN rr:=FLOAT((400-c)*900 DIV 400)*K ELSE rr:=FLOAT((c-400)*900 DIV 624)*K END;
        IF c<200 THEN gg:=FLOAT(c*900 DIV 200)*K ELSIF c<640 THEN gg:=FLOAT((640-c)*900 DIV 440)*K
        ELSE gg:=FLOAT((c-640)*900 DIV 372)*K END;
        IF c<320 THEN bb:=0.0 ELSIF c<640 THEN bb:=FLOAT((c-320)*900 DIV 320)*K ELSE bb:=900.0*K END;

        a:=rr*0.3 + gg*0.5 + bb*0.2;                 (* luminance *)     
        a:=a+0.5*(ABS(rr-a)+ABS(gg-a)+ABS(bb-a));    (* lower luminance on saturated colours *)
        a:=LUM/a;
        rr:=rr*a;
        gg:=gg*a;
        bb:=bb*a; 
        IF rr>FLOAT(LUM) THEN rr:=FLOAT(LUM) END;
        IF gg>FLOAT(LUM) THEN gg:=FLOAT(LUM) END;
        IF bb>FLOAT(LUM) THEN bb:=FLOAT(LUM) END;
      END;
    END;

  ELSE                                               (* interpolate table *)
    FILL(ADR(lut), 0C, SIZE(lut));
    FILL(ADR(tlut), 0C, SIZE(tlut));
    a:=FLOAT(atab.len-1)/FLOAT(MAXLUT);
    FOR x:=0 TO HIGH(lut) DO
      f:=a*FLOAT(x);
      y:=TRUNC(f);
      f:=f-FLOAT(y);
      IF y+2>atab.len THEN y:=atab.len-2 END;
      lut[x].rr:=intpol(atab.c[y].r, atab.c[y+1].r, f);
--WrFixed(lut[x].rr, 0, 8);
      lut[x].gg:=intpol(atab.c[y].g, atab.c[y+1].g, f);
      lut[x].bb:=intpol(atab.c[y].b, atab.c[y+1].b, f);
    END;
    FOR x:=0 TO HIGH(tlut) DO
      c:=0;
      INC(c, atab.tag[x].r); tlut[x].rr:=(LUM/255.0)*FLOAT(atab.tag[x].r);
      INC(c, atab.tag[x].g); tlut[x].gg:=(LUM/255.0)*FLOAT(atab.tag[x].g);
      INC(c, atab.tag[x].b); tlut[x].bb:=(LUM/255.0)*FLOAT(atab.tag[x].b);
      tlut[x].on:=c>0;                          (* tag is switched on *)
    END;
  END;

  amul:=amax-amin;
  IF amul<1.0 THEN amul:=1.0 END;
  amul:=1.0/amul;

  hs:=altcolstep DIV 2;
  FOR y:=0 TO HIGH(img^[0]) DO
    FOR x:=0 TO HIGH(img^) DO
      WITH img^[x][y] DO
        IF NOT ODD(b) THEN                                        (* not heaven *)
          a:=VAL(REAL, r+SHADOWLUM);                              (* luminance *)
          IF (b>0) & tlut[b DIV 2-1].on THEN                      (* tagged & tag color set *)
            WITH tlut[b DIV 2-1] DO
              r:=VAL(INTEGER, rr*a);
              g:=VAL(INTEGER, gg*a);
              b:=VAL(INTEGER, bb*a);
            END;
          ELSE
            h:=VAL(INTEGER, (VAL(REAL, g)-amin)*amul*FLOAT(HIGH(lut)));
            IF (altcolstep>0) & (h>=0) THEN h:=VAL(CARDINAL,h+hs) DIV altcolstep*altcolstep END;
            IF h<0 THEN h:=0 ELSIF h>HIGH(lut) THEN h:=HIGH(lut) END; 
            WITH lut[h] DO 
              r:=VAL(INTEGER, rr*a);
              g:=VAL(INTEGER, gg*a);
              b:=VAL(INTEGER, bb*a);
            END;
          END;
        ELSE                                                     (* heaven *)
          r:=ALTICOLORANG*LUM DIV 3;
          g:=ALTICOLORANG*LUM DIV 2;
          b:=ALTICOLORANG*LUM*3 DIV 4;
        END;
      END;
    END;
  END;
END colourize;

PROCEDURE testimage(img:pIMAGE; VAR min, max:CARDINAL);
CONST W0=2;
      V=20;
VAR x, y, w, h:CARDINAL;
BEGIN
  w:=W0+HIGH(img^) DIV 200;
  h:=w*17;
  IF (HIGH(img^[0])<=h+V) OR (HIGH(img^)<=w) THEN RETURN END;

  FOR y:=0 TO h-1 DO
    FOR x:=HIGH(img^)-(w-1) TO HIGH(img^) DO
      WITH img^[x][y+V] DO
        r:=min+(max-min)*(y DIV w) DIV 16;
        g:=r;
        b:=r;
      END;
    END;
  END;
END testimage;

PROCEDURE drawimage;
VAR panowin:PANOWIN;
--sym:SYMBOL;
BEGIN
--readsymbol(sym, "/tmp/s1.png");
  panowin.image:=image;
  panowin.flatscreen:=flatscreen;
  panowin.eye:=posa;
  panowin.horizon:=posb;
  panowin.angle:=angle;
  panowin.elevation:=elev;
  panowin.eyealt:=camalt;
  panowin.eyeog:=camog;
  panowin.yzoom:=yzoom;
  panowin.refract:=refraction;
  panowin.sunazi:=sunaz*RAD;
  panowin.sunele:=sunel*RAD;
  poirect0.lat:=PI;
  poirect0.lat:=PI;
  poirect1.long:=-PI;
  poirect1.long:=-PI;

--  panowin.azim:=azimuthl(panowin.eye, panowin.horizon)*RAD;
  panowin.hseg:=angle*RAD*0.5;
  panowin.dist:=distancel(panowin.eye, panowin.horizon);
  panowin.maxlum:=2000;
  IF imagefn[0]<>0C THEN Panorama1(panowin) END;

  IF alticolour THEN
    IF abscolmin<>0.0 THEN alticolmin:=abscolmin END; 
    IF abscolmax<>0.0 THEN alticolmax:=abscolmax END;
    colourize(panowin.image, alticolmin, alticolmax, altcoltab);
  END;

  histogram(panowin.image, alticolour, panowin.minlum, panowin.maxlum);
  makegammatab(panowin.minlum, panowin.maxlum);
  IF fonttypfoot>0 THEN
    testimage(panowin.image, panowin.minlum, panowin.maxlum);
    IF panowin.flatscreen<>eCAMERA THEN scale(panowin) END;
    WrLatLong(panowin);
  END; 
  IF trackfiles<>NIL THEN rdtracks(panowin) END;
  IF poifiles<>NIL THEN rdmountains(panowin) END;
END drawimage;

PROCEDURE camaltitude(posa, posb:POSITIONL; fromnn:BOOLEAN; VAR alta:REAL):REAL;
VAR nn, resoltx, cam:REAL;
    atwa:CARD8;
BEGIN
  IF NOT posvalidl(posa) THEN Error("need camera position") END;
  nn:=getsrtmlong(posa.lat, posa.long, 0, TRUE, resoltx, atwa, NIL);
  IF defaultocean & (nn>NOALT) THEN nn:=0.0 END;
  IF nn>=NOALT THEN closesrtmfile; Error("no altitude data at camera position") END;
  IF NOT fromnn THEN cam:=alta+nn ELSE
    cam:=alta;
    alta:=alta-nn;
  END;
  IF cam<nn THEN Error("altitude of camera below ground") END;
  nn:=getsrtmlong(posb.lat, posb.long, 0, TRUE, resoltx, atwa, NIL);
  IF defaultocean & (nn>NOALT) THEN nn:=0.0 END;
  IF nn>=NOALT THEN closesrtmfile; Error("no altitude data at view point") END;
  RETURN cam
END camaltitude;


BEGIN
  imagefn:="";
  poifiles:=NIL;
  trackfiles:=NIL;
  csvfn:="";
  csvfd:=-1;
  image:=NIL;
  xsize:=600;
  ysize:=400;
  posinval(posa);
  posinval(posb);
  camog:=10.0;
  refraction:=0.13;
  igamma:=1.0;
  srtmmaxmem:=100000000;
  heavencol.r:=cHEAVENr;
  heavencol.g:=cHEAVENg;
  heavencol.b:=cHEAVENb;
  watercol.r:=-1.0;
  watercol.g:=0.0;
  watercol.b:=0.0;
  waterglance:=0.5;
  watertransparency:=1.0;
  surfacedust:=700.0;
  elev:=0.0;
  angle:=45.0;
  yzoom:=1.0;
  flatscreen:=eCURVED;
  poisize:=5;
  minpoialt:=1.0;
  fonttyp:=0;
  fonttypfoot:=0;
  verb:=FALSE;
  poilabel:=0;
--  bicubic:=FALSE;
  rastercsv:=0;
  sunaz:=0.0;
  sunel:=-100.0;
  snowalt:=5000.0;
  slowaltslant:=500.0;
  flatwater:=FALSE;
  waveamp:=0.006;
  dusthighpass:=0.001;
  maxmountain:=8900.0;
  wavelength:=1.0;
  waterslant:=0.0;
  treesize:=0.0;
  urbanoff:=1.0;
  botanicslant:=0.0;
  ignorblack:=-0.125;
  ignorwhite:=1.0;
  whitelim:=500;
  blacklim:=500;
  desertcol.r:=-1.0;
  videolat:=1000.0;
  logfilm:=-1.0;
  newsnow:=0.0;
  sandgranularity:=1.0;
  zerotreealt:=FLOAT(ZEROTREEALT);
  fulltreealt:=FLOAT(FULLTREEALT);
  defaultocean:=FALSE;
  alticolour:=FALSE;
  alticolmin:=10000.0;
  alticolmax:=-10000.0;
  abscolmin:=0.0;
  abscolmax:=0.0;
  altcolstep:=0;
  colfn:="";
  freeshadow:=NIL;
  newshadow:=NIL;
  dropshadows:=FALSE;
  treerastersize:=TREERASTERSIZE;
debugx:=-10.0;
debugy:=10.0;
--WrFixed(fracint(-2.1, csvfd), 3, 7); WrInt(csvfd, 10); WrStrLn(" fi");
 
  Parms;
  treeraster:=40000000.0/PI2/treerastersize;
  treescramblat:=treeraster;
  treescramblong:=treeraster;

  IF xsize<20 THEN Error("xsize too less") END;
  IF ysize<20 THEN Error("ysize too less") END;
  IF (angle>=180.0) & (flatscreen<>eCURVED) THEN Error("sight > 180 deggrees needs courved screen on") END;

  readcolourfile(altcoltab, colfn);
  camalt:=camaltitude(posa, posb, videolat<90.0, camog);
  Gencrctab;
  makenoise;
--  IF imagefn[0]=0C THEN Error("need inage filename") END;

  scalesize:=fontsizex(fonttyp)*7;
--WrInt(scalesize, 5); WrStrLn("scale");
  NEW(image, xsize, ysize);
  IF image=NIL THEN Error("out of memory") END;

  drawimage;
  closesrtmfile;  
  IF csvfd>=0 THEN Close(csvfd) END;  
  IF (imagefn[0]<>0C) & (wrpng()<0) THEN Werr("png write error!"+LF) END;
  IF verb THEN WrStrLn(""); WrStr("srtmaccesses:"); WrInt(accesses, 1);
               WrStr("/"); WrInt(accessescached, 1); WrStrLn("") END;

--WrStrLn(""); WrFixed(fmod(-5.0, PI), 4, 9); WrFixed(debugx/RAD, 2, 9);WrFixed(debugy/RAD, 2, 9); WrStrLn(" xy");
--WrInt(debugxi, 12);WrInt(debugyi, 12); WrStrLn("=pixc");
END panorama.