<*+M2EXTENSIONS *>
<*-CHECKDIV *>
<*-CHECKRANGE *>
<*-COVERFLOW *>
<*-IOVERFLOW*>

<*+NOPTRALIAS*>
<*CPU="PENTIUM"*>
<*-DOREORDER *>
<*-CHECKNIL *>
<*-CHECKSET*>
<*-PROCINLINE*>
<* IF __GEN_C__ THEN *>
  <*-GENCTYPES*>
  <*+COMMENT*>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-GENDATE*>
  <*-LINENO*>
  <*-CHECKINDEX*>
  <*+GENCDIV*>
  <*-GENKRC*>
  <*-GENSIZE*>
  <*-ASSERT*>
<* ELSE *>
  <*-CHECKINDEX*>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-LINENO*>
<* END *>

IMPLEMENTATION MODULE gpspos; (* decode weather sonde gps data to position by oe5dxl *)

FROM osi IMPORT WrLn, WrStr, WrStrLn, WrInt, WrHex, Close, RdBin, WrBin, WrFixed,
                OpenRead, time; 
FROM SYSTEM IMPORT FILL, ADR, INT8, INT16, CARD16, CAST, SHIFT, CARD8, ADDRESS;
FROM math IMPORT cos,sin,atan,tan,pow,sqrt;
FROM sem IMPORT (*SEM_structAlmanac,*) SEM_ReadAlmanacDataFromFile;
FROM yuma IMPORT (*YUMA_structAlmanac,*) YUMA_ReadAlmanacDataFromFile;
FROM rinex IMPORT  RINEX_DecodeGPSNavigationFile;
FROM gps IMPORT GPS_ComputeSatellitePositionVelocityAzimuthElevationDoppler_BasedOnAlmanacData,
                (*GPS_structEphemeris,*)
                GPS_ComputeSatellitePositionVelocityAzimuthElevationDoppler_BasedOnEphmerisData;
FROM navigation IMPORT NAVIGATION_PerformClosedFormPositionSolution_FromPseuodrangeMeasurements;
FROM geodesy IMPORT GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates,
                    GEODESY_enumReferenceEllipse,
                    GEODESY_ComputeAzimuthAndElevationAnglesBetweenToPointsInTheEarthFixedFrame;
--FROM TimeConv IMPORT time;
FROM aprsstr IMPORT DateToStr, TimeToStr;

CONST
      LIGHT=299792458.0;
      FREQ=1575420000.0;
      PI=3.1415926535897932384626433832795;
      RAD=PI/180;
      PRM=LIGHT/1023000.0/1024.0;
      TESTDIL=50.0;                              (* shift sat to test influence to position *)
      WAVLEN=LIGHT/FREQ;
      MAXSPEED=500;
      SECONDSINWEEK=3600*24*7;
     
TYPE SET16=SET OF [0..15];
     SET32=SET OF [0..31];
     SATPOS=RECORD
       ok:BOOLEAN;
       x,y,z,clk,azimuth, elevation, ionocorr, range, doppler:LONGREAL;
       hdil, vdil:LONGREAL;
     END;
     SATPOSES=ARRAY[0..31] OF SATPOS;

     RESULTS=RECORD
       lat, long, heig:LONGREAL;
       hd, vd, dlat, dlong, dheig:ARRAY[0..3] OF LONGREAL;        (* h and v dilatation per sat *)
       qsumv, qsumh:LONGREAL;
       satset:SET16;
       res:SET16;
     END;

    COMMONALMANACH=RECORD
      treal,
      toa:CARDINAL;          (* almanac time of applicability (reference time [s]*)
      week,                  (* 10 bit gps week 0-1023 (user must account for week rollover) [week]*)
      prn,                   (* GPS prn number *)
      reserved,              (* reserved *)
      svn:CARD16;            (* Satellite vehicle number *)
      ura,                   (* User Range Accuracy lookup code, [0-15], see p. 83 GPSICD200C, 0 is*)
      health,                (* 0=healthy, unhealthy otherwise  [], subframe 4 and 5, page 25 six-b*)
      config_code:CARD8;      (* configuration code   [], if >=9 Anti-Spoofing is on *)

  (* this inicator is not part of the SEM standard but is added by the user if known *)
      is_af0_af1_high_precision:CHAR;  (* indicates precision of af0 and af1 [1=high precision,0=lo*)

      ecc,          (* eccentricity *)
      i0,           (* orbital inclination at reference time                        [rad] *)
      omegadot,     (* rate of right ascension                                      [rad/s] *)
      sqrta,        (* square root of the semi-major axis                           [m^(1/2)]*)
      omega0,       (* longitude of ascending node of orbit plane at weekly epoch   [rad] *)
      w,            (* argument of perigee                                          [rad] *)
      m0,           (* mean anomaly at reference time                               [rad] *)
      af0,          (* polynomial clock correction coefficient (clock bias)         [s],   Note: pa*)
      af1,
      af2:LONGREAL; (* polynomial clock correction coefficient (clock drift)        [s/s], Note: pa*)


      toe,                (* reference time ephemeris (0-604800) *)
      toc:CARDINAL;       (* reference time (clock)   (0-604800) *)
      iodc:CARD16;        (* 10 bit issue of data (clock) *)
      iode,               (* 8 bit  issue of data (ephemeris) *)
      alert_flag,         (* 1 = URA may be worse than indicated *)
      anti_spoof,         (* anti-spoof flag from 0=off, 1=on *)
      code_on_L2,         (* 0=reserved, 1=P code on L2, 2=C/A on L2 *)
      L2_P_data_flag,     (* flag indicating if P is on L2 1=true *)
      fit_interval_flag:CHAR;  (* fit interval flag (four hour interval or longer) 0=4 fours *)
      age_of_data_offset, (* age of data offset *)
      tow_week:CARD16;    (* The week corresponding to tow (0-1024+). Can be one week l *)
      tow:CARDINAL;       (* The time of week derived formt the Z-count in the Hand Ove *)

(* clock parameters *)
      tgd:LONGREAL;       (* group delay *)


(* ephemeris parameters *)
      delta_n,    (* mean motion difference from computed value *)
      idot,       (* rate of inclination angle *)
      cuc,        (* amplitude of the cosine harmonic correction term to the argument o *)
      cus,        (* amplitude of the sine harmonic correction term to the argument of *)
      crc,        (* amplitude of the cosine harmonic correction term to the orbit radi *)
      crs,        (* amplitude of the sine harmonic correction term to the orbit radius *)
      cic,        (* amplitude of the cosine harmonic correction term to the angle of i  *)
      cis:LONGREAL; (* amplitude of the sine harmonic correction term to the angle of inc *)

    END;

  structEphemeris=RECORD
    toe,                (* reference time ephemeris (0-604800) *)
    toc:CARDINAL;       (* reference time (clock)   (0-604800) *)
    prn,                (* GPS PRN number *)
    week,               (* 10 bit gps week 0-1023 (user must account for week rollove *)
    iodc,               (* 10 bit issue of data (clock) *)
    reserved1:CARD16;   (* reserved bytes *)
    iode,               (* 8 bit  issue of data (ephemeris) *)
    health,             (* 6 bit health parameter, 0 if healthy, unhealth othersize *)
    alert_flag,         (* 1 = URA may be worse than indicated *)
    anti_spoof,         (* anti-spoof flag from 0=off, 1=on *)
    code_on_L2,         (* 0=reserved, 1=P code on L2, 2=C/A on L2 *)
    L2_P_data_flag,     (* flag indicating if P is on L2 1=true *)
    fit_interval_flag,  (* fit interval flag (four hour interval or longer) 0=4 fours *)
    ura:CHAR;           (* User Range Accuracy lookup code, 0 is excellent, 15 is use *)
    age_of_data_offset, (* age of data offset *)
    tow_week:CARD16;    (* The week corresponding to tow (0-1024+). Can be one week l *)
    tow:CARDINAL;       (* The time of week derived formt the Z-count in the Hand Ove *)

(* clock parameters *)
    tgd,          (* group delay *)
    af2,          (* polynomial clock correction coefficient (rate of clock drift) *)
    af1,          (* polynomial clock correction coefficient (clock drift) *)
    af0:LONGREAL; (* polynomial clock correction coefficient (clock bias) *)

(* ephemeris parameters *)
    m0,         (* mean anomaly at reference time *)
    delta_n,    (* mean motion difference from computed value *)
    ecc,        (* eccentricity *)
    sqrta,      (* square root of the semi-major axis *)
    omega0,     (* longitude of ascending node of orbit plane at weekly epoch *)
    i0,         (* inclination angle at reference time *)
    w,          (* argument of perigee *)
    omegadot,   (* rate of right ascension *)
    idot,       (* rate of inclination angle *)
    cuc,        (* amplitude of the cosine harmonic correction term to the argument o *)
    cus,        (* amplitude of the sine harmonic correction term to the argument of *)
    crc,        (* amplitude of the cosine harmonic correction term to the orbit radi *)
    crs,        (* amplitude of the sine harmonic correction term to the orbit radius *)
    cic,        (* amplitude of the cosine harmonic correction term to the angle of i  *)
    cis:LONGREAL; (* amplitude of the sine harmonic correction term to the angle of inc *)
  END;

  structKlobuchar=RECORD
    isValid:CARD16; (* Is this structure valid for use 1=YES, 0=NO.*)
    week:CARD16;    (* The GPS week corresponding to the correction parameters [weeks]. *)
    tow:CARDINAL;   (* The GPS time of week corresponding to the correction parameters [s]. *)
    alpha0,         (* coefficients of a cubic equation representing the amplitude of the ve *)
    alpha1,         (* coefficients of a cubic equation representing the amplitude of the ve *)
    alpha2,         (* coefficients of a cubic equation representing the amplitude of the ve *)
    alpha3,         (* coefficients of a cubic equation representing the amplitude of the ve *)
    beta0,          (* coefficients of a cubic equation representing the period of the model *)
    beta1,          (* coefficients of a cubic equation representing the period of the model *)
    beta2,          (* coefficients of a cubic equation representing the period of the model *)
    beta3:LONGREAL; (* coefficients of a cubic equation representing the period of the model *)
  END;


VAR rinexklobuchar:structKlobuchar; 
    semok, yumaok, rinexok:BOOLEAN;
    stats:ARRAY[0..499] OF RESULTS;
    rangcor:ARRAY[0..31] OF INTEGER;
    calm:ARRAY[0..31] OF COMMONALMANACH;

PROCEDURE sqr(x:LONGREAL):LONGREAL;
BEGIN RETURN x*x END sqr;

PROCEDURE degtostr(d:REAL; lat:BOOLEAN; form:CHAR; VAR s:ARRAY OF CHAR);
CONST Z=ORD("0");
VAR n,i:CARDINAL;
BEGIN
  IF HIGH(s)<11 THEN s[0]:=0C; RETURN END;

  IF form="2" THEN i:=7 ELSIF form="3" THEN i:=8 ELSE i:=9 END;
  IF d<0.0 THEN
    d:=-d;
    IF lat THEN s[i]:="S" ELSE s[i+1]:="W" END;
  ELSIF lat THEN s[i]:="N" ELSE s[i+1]:="E" END;

  IF form="2" THEN                                          (* DDMM.MMNDDMM.MME *)
    n:=TRUNC(d*(6000*180/PI)+0.5);
    s[0]:=CHR(n DIV 600000 MOD 10+Z);
    i:=ORD(NOT lat);
    s[i]:=CHR(n DIV 60000 MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV 6000 MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV 1000 MOD 6+Z); INC(i);
    s[i]:=CHR(n DIV 100 MOD 10+Z); INC(i);
    s[i]:="."; INC(i);
    s[i]:=CHR(n DIV 10 MOD 10+Z); INC(i);
    s[i]:=CHR(n MOD 10+Z); INC(i);

  ELSIF form="3" THEN                                        (* DDMM.MMMNDDMM.MMME *)
    n:=TRUNC(d*(60000*180/PI)+0.5);
    s[0]:=CHR(n DIV 6000000 MOD 10+Z);
    i:=ORD(NOT lat);
    s[i]:=CHR(n DIV 600000 MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV 60000 MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV 10000 MOD 6+Z); INC(i);
    s[i]:=CHR(n DIV 1000 MOD 10+Z); INC(i);
    s[i]:="."; INC(i);
    s[i]:=CHR(n DIV 100 MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV 10 MOD 10+Z); INC(i);
    s[i]:=CHR(n MOD 10+Z); INC(i);

  ELSE                                                       (* DDMMSS *)
    n:=TRUNC(d*(60*60*180/PI)+0.5);
    s[0]:=CHR(n DIV (60*6000) MOD 10+Z);
    i:=ORD(NOT lat);
    s[i]:=CHR(n DIV (60*600) MOD 10+Z); INC(i);
    s[i]:=CHR(n DIV (60*60) MOD 10+Z); INC(i);
    s[i]:="d"; INC(i);
    s[i]:=CHR(n DIV 600 MOD 6+Z); INC(i);
    s[i]:=CHR(n DIV 60 MOD 10+Z); INC(i);
    s[i]:="'"; INC(i);
    s[i]:=CHR(n DIV 10 MOD 6+Z); INC(i);
    s[i]:=CHR(n MOD 10+Z); INC(i);
    s[i]:='"'; INC(i);
  END;

  INC(i);
  s[i]:=0C;
END degtostr;


PROCEDURE wgs84(lat, long, heig:LONGREAL; VAR x,y,z:LONGREAL);       (* wgs84 ecef *)
CONST
--      A=6378217; 
      A=6378137;
--      B=6356752;
--      F=(A-B)/A;
--      E=2*F-F*F;
      E=0.081819190842522; 

VAR n,sl,h:LONGREAL;
BEGIN
  sl:=sin(lat);
  n:=A/sqrt(1.0-E*E*sl*sl);
  h:=heig+n;
  z:=(n*(1.0-E*E)+heig)*sl;
  y:=h*sin(long)*cos(lat);
  x:=h*cos(long)*cos(lat);
END wgs84;


PROCEDURE get4sats(sats-:SATPOSES; satnum-:ARRAY OF CARDINAL;
                   dil:CARDINAL; dilm:LONGREAL; VAR lat, long, heig:LONGREAL):INTEGER;
CONST K=0.0;
VAR i, sat:CARDINAL;
    clk, x,y,z:ARRAY[0..3] OF LONGREAL;
    rx_clock_bias (*, clock_drift, satVx,satVy,satVz,azimuth,elevation,doppler*):LONGREAL;
    ret, prn:INTEGER;
    chkmask:SET32;
    dils:ARRAY[0..3] OF LONGREAL;

BEGIN
  chkmask:=SET32{};
  FOR i:=0 TO HIGH(dils) DO
    dils[i]:=0.0;
    prn:=satnum[i]-1;
    IF (prn<0) OR (prn>31) OR (prn IN chkmask) THEN RETURN -1 END;      (* two times same sat *)

    INCL(chkmask, prn);
  END;
  IF dil>0 THEN dils[dil-1]:=dilm END;  
  ret:=NAVIGATION_PerformClosedFormPositionSolution_FromPseuodrangeMeasurements(
   sats[satnum[0]].range+dils[0](*+K*sats[satnum[0]].drydelay*),
   sats[satnum[1]].range+dils[1](*+K*sats[satnum[1]].drydelay*),
   sats[satnum[2]].range+dils[2](*+K*sats[satnum[2]].drydelay*),
   sats[satnum[3]].range+dils[3](*+K*sats[satnum[3]].drydelay*),
   sats[satnum[0]].clk, sats[satnum[1]].clk, sats[satnum[2]].clk, sats[satnum[3]].clk,
   sats[satnum[0]].x, sats[satnum[1]].x, sats[satnum[2]].x, sats[satnum[3]].x,
   sats[satnum[0]].y, sats[satnum[1]].y, sats[satnum[2]].y, sats[satnum[3]].y,
   sats[satnum[0]].z, sats[satnum[1]].z, sats[satnum[2]].z, sats[satnum[3]].z,
   lat, long, heig, rx_clock_bias);
--WrFixed(rx_clock_bias, 2,12); WrStrLn("=bias");
   IF ret>500 THEN ret:=500 END;
   IF heig<-10000.0 THEN ret:=0 END; 
   RETURN ret

<* IF NOT __GEN_C__ THEN *>

EXCEPT
  RETURN -1
<* END *>
END get4sats;

PROCEDURE qdist(dlat, dlong:LONGREAL):LONGREAL; (* horizontal distance square rad*)
VAR c:LONGREAL;
BEGIN
  c:=cos(dlat)*dlong;
  RETURN dlat*dlat + c*c
END qdist;


PROCEDURE satposit(VAR sat:SATPOS; satnum:CARDINAL; myx, myy, myz:LONGREAL;
                   userweek, weekms:CARDINAL);
VAR 
    rx_clock_bias, clock_drift, clock_correction, satVx,satVy,satVz:LONGREAL;

BEGIN
  sat.ok:=FALSE;

  IF calm[satnum].prn>0 THEN
    WITH calm[satnum] DO
      IF rinexok THEN
        GPS_ComputeSatellitePositionVelocityAzimuthElevationDoppler_BasedOnEphmerisData(
        myx, myy, myz, userweek, LFLOAT(weekms)*0.001,
        tow_week MOD 1024,
        toe, toc, af0, af1, af2, tgd, m0, delta_n, ecc, sqrta, omega0, i0, w, omegadot, idot,
        cuc, cus, crc, crs, cic, cis,
        sat.clk, clock_drift, sat.x, sat.y, sat.z, satVx, satVy, satVz,
        sat.azimuth, sat.elevation, sat.doppler);

--WrInt(ORD(fit_interval_flag), 4); WrInt(tow, 10); WrInt(weekms DIV 1000, 10); WrStrLn("=fit tow tsat");
 

      ELSE
        GPS_ComputeSatellitePositionVelocityAzimuthElevationDoppler_BasedOnAlmanacData(
        myx, myy, myz, userweek, LFLOAT(weekms)*0.001, FLOAT(toa),
        week, prn, ecc, i0, omegadot, sqrta, omega0, w, m0, af0, af1,
        sat.clk, clock_drift, sat.x, sat.y, sat.z, satVx,satVy,satVz,
        sat.azimuth, sat.elevation, sat.doppler);
      END;
--WrInt(toe, 9); WrInt((time()-7200+3600*24*3) MOD 604800, 9); WrInt(weekms DIV 1000, 9); WrStrLn("=toe");
(*
WrInt(prn, 2);
WrFixed(ecc, 5, 9);
WrFixed(m0, 5, 9);
WrFixed(sqrta, 5, 9);
WrFixed(omega0, 5, 9);
WrFixed(i0, 5, 9);
WrFixed(w, 5, 9);
WrFixed(omegadot, 5, 9);

--WrInt(toe, 5);
--WrInt(toc, 5);
WrInt(tow_week, 5);
WrFixed(sat.azimuth/RAD, 1, 7);
WrFixed(sat.elevation/RAD, 1, 6);
WrStrLn("=p a e"); 
*)
    END;    
    sat.ok:=TRUE;
  END;

--WrInt(alm[satnum].prn, 3); WrFixed(clock_drift, 7, 10); WrStrLn("=clock_drift");


END satposit;


PROCEDURE lowele(e:LONGREAL):LONGREAL;
BEGIN
  IF e<0.01 THEN e:=0.01 END;
  RETURN 1.0/sin(e)
END lowele; 


PROCEDURE stat4sats(VAR sats:SATPOSES; satnums-:ARRAY OF CARDINAL;
                    try:CARDINAL):BOOLEAN;
VAR i:CARDINAL;
    ret:INTEGER;
    bx,by,bz, dil, ele, azi: LONGREAL;

BEGIN
  WITH stats[try] DO
    IF get4sats(sats, satnums, 0, 0.0, lat, long, heig)>0 THEN

      wgs84(lat, long, heig, bx, by, bz);                      (* for elevation *)
      FOR i:=0 TO 3 DO

        dil:=TESTDIL;
        IF GEODESY_ComputeAzimuthAndElevationAnglesBetweenToPointsInTheEarthFixedFrame(
           GEODESY_REFERENCE_ELLIPSE_WGS84, bx, by, bz,
           sats[satnums[i]].x, sats[satnums[i]].y, sats[satnums[i]].z, ele, azi)
        THEN dil:=dil+ lowele(ele) * (TESTDIL*0.1) END;        (* make bad hdop on low elevation *)

        IF get4sats(sats, satnums, i+1, dil, dlat[i], dlong[i], dheig[i])>0 THEN
          hd[i]:=sqrt(qdist(lat-dlat[i], long-dlong[i]))*(20000000/PI);
          vd[i]:=ABS(heig-dheig[i]);
          sats[satnums[i]].hdil:=sats[satnums[i]].hdil+hd[i];
          sats[satnums[i]].vdil:=sats[satnums[i]].vdil+vd[i];

        ELSE RETURN FALSE END;

      END;
      RETURN TRUE

    END;
  END;
  RETURN FALSE
END stat4sats;


PROCEDURE satposits(VAR satsp:SATPOSES; myx, myy, myz:LONGREAL; userweek, weekms:CARDINAL);
VAR i:CARDINAL;
    rx_clock_bias, clock_drift, satVx,satVy,satVz:LONGREAL;

BEGIN
  FILL(ADR(satsp), 0C, SIZE(satsp));
  FOR i:=0 TO HIGH(satsp) DO
    satposit(satsp[i], i, myx, myy, myz, userweek, weekms);
  END;
END satposits;

(*
PROCEDURE tropomodel(VAR satpos:SATPOS; userlat, userhigh:LONGREAL; time:CARDINAL);
VAR    zenith_dry_delay, zenith_wet_delay :LONGREAL;
       dayofyear:CARDINAL;
BEGIN
  WITH satpos DO
    IF elevation>0.0 THEN
      dayofyear:=time DIV (3600*24 DIV 4) MOD (365*4+1) DIV 4;
      TROPOSPHERE_DetermineZenithDelayValues_WAAS_Model(userlat, userhigh, dayofyear,
      zenith_dry_delay, zenith_wet_delay);

      TROPOSPHERE_GetDryAndWetDelay_UsingThe_UNBabc_MappingFunction(
      zenith_dry_delay, zenith_wet_delay, elevation, userlat, userhigh,
      drydelay, wetdelay);

--WrFixed(elevation/RAD, 1,6); WrFixed(wetdelay, 5,12); WrFixed(drydelay, 5,12);
--WrFixed(zenith_dry_delay, 5,12); WrFixed(zenith_wet_delay, 5,12); WrStrLn("");
    END;
  END;

END tropomodel;
*)

PROCEDURE killexo(VAR stats:ARRAY OF RESULTS; probs:CARDINAL; VAR devhsum, devvsum:REAL);
CONST MAXHMUL=2.0;
      MAXVMUL=4.0;
VAR i:CARDINAL;
    c, im:CARDINAL;
    dx, dy:LONGREAL;
    dev, max, cosy:REAL;
    ok:BOOLEAN;
BEGIN
  REPEAT
    ok:=TRUE;
    c:=0;
    dx:=0.0;
    dy:=0.0; 
    FOR i:=0 TO probs-1 DO                         (* median *)
      IF stats[i].res*SET16{0,3}=SET16{} THEN
        dx:=dx+stats[i].long;
        dy:=dy+stats[i].lat;
        INC(c);
      END;
    END;
    devhsum:=0.0;
    IF c>1 THEN                                    (* filter min 3 measurements only *) 
      dx:=dx/LFLOAT(c);
      dy:=dy/LFLOAT(c);
      cosy:=cos(VAL(REAL,dy));
      max:=-1.0;
      FOR i:=0 TO probs-1 DO                       (* deviation and find badest *)
        IF stats[i].res*SET16{0,3}=SET16{} THEN
          dev:=sqr(VAL(REAL, stats[i].long-dx)*cosy) + sqr(VAL(REAL, stats[i].lat-dy));   (* error distance ^2 *)
--WrFixed(sqrt(dev)/RAD/360.0*40000000.0, 1,6); WrStrLn("=dev");
          devhsum:=devhsum + dev;  
          IF dev>max THEN max:=dev; im:=i END;
        END;
      END;
      IF (max>=0.0) & (c>2) THEN INCL(stats[im].res, 0) END;      (* mark as bad *)
      devhsum:=devhsum/FLOAT(c);                   (* median deviation rad^2*)
      IF max>devhsum*MAXHMUL THEN ok:=FALSE END;   (* do until all in limit or 1 remains *)
    END;
  UNTIL ok;

  IF devhsum>0.0 THEN devhsum:=sqrt(devhsum)*(40000000.0/RAD/360.0) END;   (* meter *)    

  REPEAT
    ok:=TRUE;
    c:=0;
    dx:=0.0;
    FOR i:=0 TO probs-1 DO
      IF stats[i].res*SET16{2,4}=SET16{} THEN
        dx:=dx+stats[i].heig;
        INC(c);
      END;
    END;

    devvsum:=0.0;
    IF c>1 THEN
      dx:=dx/LFLOAT(c);
      max:=-1.0;
      FOR i:=0 TO probs-1 DO
        IF stats[i].res*SET16{2,4}=SET16{} THEN
          dev:=sqr(VAL(REAL, stats[i].heig-dx));
          devvsum:=devvsum + dev;
          IF dev>max THEN max:=dev; im:=i END;
        END;
      END;
      IF (max>=0.0) & (c>2) THEN INCL(stats[im].res, 2) END;
      devvsum:=devvsum/FLOAT(c);                   (* median alt deviation m^2 *)
      IF max>devvsum*MAXVMUL THEN ok:=FALSE END;
    END;
  UNTIL ok;

  IF devvsum>0.0 THEN devvsum:=sqrt(devvsum) END;
  
END killexo;

PROCEDURE killdop(VAR stats:ARRAY OF RESULTS; probs:CARDINAL);
VAR i:CARDINAL;
BEGIN
  FOR i:=0 TO probs-1 DO
    IF stats[i].qsumh>100000.0 THEN INCL(stats[i].res, 3) END;
    IF stats[i].qsumv>400000.0 THEN INCL(stats[i].res, 4) END;
  END;
END killdop;


PROCEDURE median(stats-:ARRAY OF RESULTS; tries:CARDINAL; VAR lat, long, heig:LONGREAL; VAR cnt:CARDINAL);
VAR i, cz:CARDINAL;
BEGIN
  lat:=0.0;
  long:=0.0;
  heig:=0.0;
  cnt:=0;
  cz:=0;
  FOR i:=0 TO tries-1 DO
    IF stats[i].res*SET16{0,3}=SET16{} THEN
      long:=long+stats[i].long;
      lat :=lat +stats[i].lat;
      INC(cnt);
    END;

    IF stats[i].res*SET16{2,4}=SET16{} THEN
      heig:=heig+stats[i].heig;
      INC(cz);
    END;
  END;
  IF cnt>0 THEN
    lat :=lat /LFLOAT(cnt);
    long:=long/LFLOAT(cnt);
  ELSE
    lat :=0.0;
    long:=0.0;
  END;
  IF cz>0 THEN heig:=heig/LFLOAT(cz) ELSE heig:=0.0 END;
END median;


PROCEDURE showstats(stats-:ARRAY OF RESULTS; probs:CARDINAL; errh, errv:REAL;
                    restcnt:CARDINAL; full:BOOLEAN);
VAR i, p:CARDINAL;
    h:ARRAY[0..30] OF CHAR;
    las, los, la, lo:LONGREAL;
BEGIN
  IF full THEN
    WrStr(" probs:");WrInt(probs, 1);WrStrLn("");
    FOR p:=0 TO probs-1 DO
      FOR i:=0 TO 11 DO WrInt(ORD(i IN stats[p].satset), 1) END;
      WrStr(" ");
      degtostr(stats[p].lat, TRUE, "2", h); WrStr(h);  WrStr(" ");
      degtostr(stats[p].long, FALSE, "2", h); WrStr(h);
      WrInt(VAL(INTEGER, stats[p].heig), 7);
      FOR i:=0 TO 3 DO 
        WrFixed(stats[p].hd[i], 0, 5);
        WrFixed(stats[p].vd[i], 0, 5);
      END;
      WrFixed(stats[p].qsumh, 0, 8);
      WrFixed(stats[p].qsumv, 0, 8); WrStr(" ");
--    WrFixed(stats[p].minele/RAD, 0, 3); WrStr(" ");
      FOR i:=0 TO 4 DO WrInt(ORD(i IN stats[p].res), 1) END;
      WrStrLn("");
    END;
  END; 
(*
  la:=0.0;
  lo:=0.0;
  i:=0;
  FOR p:=0 TO probs-1 DO la:=la+stats[p].lat; lo:=lo+stats[p].long; INC(i) END;
  la:=la/LFLOAT(i);
  lo:=lo/LFLOAT(i);
  las:=0.0;
  los:=0.0;
  FOR p:=0 TO probs-1 DO las:=las+ABS(stats[p].lat-la); los:=los+ABS(stats[p].long-lo) END;
  WrStr("s=");WrFixed(las/LFLOAT(i)*100000.0, 2, 8); WrStr(" "); WrFixed(los/LFLOAT(i)*100000.0, 2, 8); WrStrLn("");
*)
  WrStr("n="); WrInt(restcnt, 3) ;WrStr(" s=");
  WrFixed(errh, 1,0); WrFixed(errv, 1,6); WrStrLn("");
END showstats;


PROCEDURE speedrange(VAR sats:SATS; satcnt:CARDINAL; VAR satspos:SATPOSES); (* shift ranges with speed *)

  PROCEDURE median(bad:CARDINAL):LONGREAL;     (* median except from bad and badspeeds *)
  VAR i,n:CARDINAL;
      m:LONGREAL;
  BEGIN
    m:=0.0;
    n:=0;
    FOR i:=0 TO satcnt-1 DO
      IF (i<>bad) & NOT sats[i].badspeed THEN m:=m + sats[i].userspeed; INC(n) END;
    END;
    IF n>0 THEN m:=m/LFLOAT(n) END;
    RETURN m
  END median;

  PROCEDURE peak(med:LONGREAL; bad:CARDINAL):CARDINAL;  (* max deviation except from bad and badspeeds *)
  VAR i,b:CARDINAL;
      m,a:LONGREAL;
  BEGIN
    m:=0.0;
    b:=MAX(CARDINAL);
    FOR i:=0 TO satcnt-1 DO
      IF (i<>bad) & NOT sats[i].badspeed THEN
        a:=ABS(sats[i].userspeed-med);
        IF a>=m THEN m:=a; b:=i END;
      END;
    END;
    RETURN b
  END peak;

CONST MINSATS=4;
      MAXSPEED=200;
VAR i, bad, bad1, goodsats:CARDINAL;
    med, med1:LONGREAL;
BEGIN
  IF satcnt>=MINSATS THEN
    FOR i:=0 TO satcnt-1 DO
      WITH sats[i] DO
        userspeed:=LFLOAT(rang1)*(WAVLEN/256.0)-satspos[almidx].doppler;  (* doppler speed sat - user *)
        badspeed:=FALSE;
      END;
    END;

    goodsats:=satcnt;
    LOOP
--WrStr("speeds");
      med:=median(MAX(CARDINAL));          (* find sat oszillator deviation *)
      bad:=peak(med, MAX(CARDINAL));       (* find bad synchronized doppler for a sat *)
      IF bad>=satcnt THEN EXIT END;
--WrInt(bad, 3);

      med1:=median(bad);                   (* median without bad sat *)
      IF ABS(sats[bad].userspeed-med1)<=MAXSPEED THEN
        bad1:=peak(med1, bad);               (* second badest sat *)
--WrInt(bad1, 3);
        IF (bad1>=satcnt) OR (ABS(sats[bad].userspeed-med1)*0.25<ABS(sats[bad1].userspeed-med1))
        THEN EXIT END;                       (* if badest sat deviates less than k* second badest *)

      END;
      sats[bad].badspeed:=TRUE;
      med:=med1;
--WrStrLn("");

      FOR i:=0 TO satcnt-1 DO 
--WrFixed(sats[i].userspeed, 1, 7);
        sats[i].userspeed:=sats[i].userspeed-med;
--WrFixed(sats[i].userspeed, 1, 7);
--WrInt(ORD(sats[i].badspeed), 2); WrStrLn("");
      END;
      DEC(goodsats);
      IF goodsats<=MINSATS THEN EXIT END;
    END;

  END;
END speedrange;


PROCEDURE neardist(lat1, long1, lat2, long2:LONGREAL):LONGREAL;
VAR x,y:LONGREAL;
BEGIN
  x:=(long1-long2)*cos((lat1+lat2)*0.5);
  y:=lat1-lat2;
  RETURN 6370000.0*sqrt(x*x+y*y)
END neardist;

PROCEDURE azimuth(lat1, long1, lat2, long2:LONGREAL):LONGREAL;   (* degrees *)
VAR h, ldiff, sinld, cl1, cl2:LONGREAL;
BEGIN
  ldiff:=long2-long1;
  sinld:=sin(ldiff);
  cl1:=cos(lat1);
  cl2:=cos(lat2);

  IF (ldiff=0.0) OR (cl1=0.0) OR (cl2=0.0) THEN
    IF lat1<=lat2 THEN RETURN 0.0 ELSE RETURN 180.0 END;
  ELSE
    h:=atan(cl1*((tan(lat1)*cos(ldiff)-tan(lat2))/sinld))*(1.0/RAD);
    IF ldiff<0.0 THEN RETURN h+270.0 ELSE RETURN h+90.0 END;
  END;

<* IF NOT __GEN_C__ THEN *>

EXCEPT
  RETURN 0.0
<* END *>
END azimuth;

PROCEDURE dirmed(VAR dirsum:LONGREAL; dir:LONGREAL; cnt:CARDINAL); (* angle median in degrees *)
VAR d:LONGREAL;
BEGIN
--WrFixed(dir, 0,0); WrFixed(dirsum/(LFLOAT(cnt)+0.0001), 0,0);
  IF cnt=0 THEN dirsum:=dir;
  ELSE
    dirsum:=dirsum+dir;
--WrFixed(dir-dirsum/LFLOAT(cnt+1), 0,0);
    IF ABS(dir-dirsum/LFLOAT(cnt+1))>180.0 THEN dirsum:=dirsum+360.0 END;
    d:=dirsum/LFLOAT(cnt+1);
    IF d>=360.0 THEN
      d:=d-360.0;
      dirsum:=d*LFLOAT(cnt+1); 
    END;
  END;  
--WrFixed(dirsum/(LFLOAT(cnt+1)+0.0001), 0,0); WrStrLn("=dir");
END dirmed;


(*
PROCEDURE ionomodel(VAR sat:SATPOS; klobuchar:GNSS_structKlobuchar;
          mylat, mylong, myhigh:LONGREAL; tow:CARDINAL);
BEGIN
  WITH klobuchar DO
    IF NOT IONOSPHERE_GetL1KlobucharCorrection
       (alpha0, alpha1, alpha2, alpha3, beta0,  beta1,  beta2,  beta3,
       mylat, mylong, sat.elevation, sat.azimuth, LFLOAT(tow), sat.ionocorr)
    THEN sat.ionocorr:=0.0 END;
--WrFixed(sat.ionocorr, 1, 6); WrFixed(mylat/RAD, 3, 7); WrFixed(mylong/RAD, 3, 7);
--WrFixed(sat.azimuth/RAD, 1, 6);WrFixed(sat.elevation/RAD, 1, 6);WrStrLn("=iono"); 


  END;
END ionomodel;
*)


PROCEDURE PMUL(n:INTEGER):LONGREAL;
BEGIN
  RETURN VAL(LONGREAL, MAX(INTEGER)-n)*PRM;

<* IF NOT __GEN_C__ THEN *>

EXCEPT
  WrStrLn(" lfloat exception");
  RETURN 0.0
<* END *>
END PMUL;

PROCEDURE wrdate(t:CARDINAL);
VAR s:ARRAY[0..30] OF CHAR;
BEGIN
  DateToStr(t, s); WrStr(s);
END wrdate;


PROCEDURE getposit(weekms:CARDINAL; VAR systime:CARDINAL; sats:SATS; mylat, mylong,
                   myhigh:LONGREAL; VAR lat, long, heig, speed, dir, climb:LONGREAL;
                   VAR hrms, vrms:REAL; VAR goodsats:CARDINAL):INTEGER;

CONST ZEROTIME=315964800;
      WEEK=3600*24*7;
      HALFWEEK=WEEK DIV 2;

VAR myx, myy, myz, bx, by, bz, mlat, mlong, malt, min, vlat, vlong, vheig, vdx, aa, hprm:LONGREAL;
    satnums:ARRAY[0..3] OF CARDINAL;
    ret:INTEGER;
    gpstime, userweek, weeks, cnt, i, j, bit, satcnt, tries, speeds, restcnt:CARDINAL;
    comb:CARD16;
    satspos:SATPOSES;
    h:ARRAY[0..30] OF CHAR;
BEGIN
  weeks:=weekms DIV 1000 MOD WEEK;                       (* second in week from gps *)
  gpstime:=systime-ZEROTIME;                             (* absolute gps time from system clock *)
  i:=(weeks + WEEK - gpstime MOD WEEK) MOD WEEK;         (* gps is later *)
  IF i<HALFWEEK THEN INC(i, WEEK) END; 
  gpstime:=gpstime+i-WEEK;                               (* system time corrected to gps time *)
  systime:=gpstime + ZEROTIME;                           (* return corrected systime *)
  userweek:=gpstime DIV (3600*24*7) MOD 1024; 
--WrInt(calm[0].week MOD 1024, 6); WrInt(userweek, 6); WrInt(weeks, 16);WrInt(gpstime MOD WEEK, 16);WrInt(i-WEEK, 16); WrStrLn("=weeks");


  wgs84(mylat, mylong, myhigh, myx, myy, myz);

  FILL(ADR(stats), 0C, SIZE(stats));
  FILL(ADR(satspos), 0C, SIZE(satspos));
  satposits(satspos, myx, myy, myz, userweek, weekms);
--corrrang(sats);

  satcnt:=0;
  FOR i:=0 TO HIGH(sats) DO 
    IF (sats[i].prn<>0) & (sats[i].rang<MAX(INTEGER)) THEN
      j:=0;
      
      LOOP
        IF calm[j].prn=sats[i].prn THEN
          sats[satcnt]:=sats[i];
          sats[satcnt].almidx:=j;
(*
          IF rinexok THEN
            ionomodel(satspos[j], rinexklobuchar, mylat, mylong, myhigh, weeks);
          END;
*)
--        tropomodel(satspos[j], mylat, myhigh, systime);
          INC(satcnt);
          EXIT
        END;      
        INC(j);
        IF j>HIGH(calm) THEN EXIT END;   (* sat not in almanach *)
      END;
    END;
  END;

  tries:=0;

  hrms:=0.0;
  ret:=-1;
  mlat:=0.0;
  mlong:=0.0;
  malt:=0.0;
  IF satcnt>=4 THEN  
    comb:=0;
    REPEAT
      cnt:=0;
      FOR bit:=0 TO satcnt-1 DO IF (bit IN CAST(SET16, comb)) THEN INC(cnt) END; END;
      IF cnt=4 THEN 
        cnt:=0;
        FOR bit:=0 TO satcnt-1 DO
          IF (bit IN CAST(SET16, comb)) THEN 
            satspos[sats[bit].almidx].range:=PMUL(sats[bit].rang);
--            satspos[sats[bit].almidx].range:=VAL(REAL, VAL(CARDINAL,MAX(INTEGER)-sats[bit].rang))*PRM;
-- (LFLOAT(ORD(ODD(weeks))*2)-1.0)*satspos[sats[bit].almidx].ionocorr;
            satnums[cnt]:=sats[bit].almidx; 
            INC(cnt);
          END;
        END; 
        IF stat4sats(satspos, satnums, tries) THEN
          stats[tries].qsumh:=0.0;
          stats[tries].qsumv:=0.0;
          FOR i:=0 TO 3 DO 
            stats[tries].qsumh:=stats[tries].qsumh+sqr(stats[tries].hd[i]);
            stats[tries].qsumv:=stats[tries].qsumv+sqr(stats[tries].vd[i]);
            IF (stats[tries].hd[i]=0.0) OR (stats[tries].vd[i]=0.0)
            THEN stats[tries].qsumh:=0.0; stats[tries].qsumv:=0.0 END;
          END;
          stats[tries].satset:=CAST(SET16, comb);
          INC(tries);
        END;

      END;
      INC(comb);
    UNTIL (satcnt IN CAST(SET16, comb));

  END; 

WrStr(" "); wrdate(systime);
WrInt(satcnt, 3); WrStr("=sats"); WrInt(tries, 5); WrStr("=tries ");

  IF tries>0 THEN
    killdop(stats, tries);
    killexo(stats, tries, hrms, vrms);
    median(stats, tries, lat, long, heig, restcnt);   (* final position *) 
    showstats(stats, tries, hrms, vrms, restcnt, FALSE);
    ret:=0;
  END;


  goodsats:=0;
  IF ret>=0 THEN
-----speed
    IF (lat=0.0) & (long=0.0) THEN wgs84(mylat, mylong, myhigh, bx, by, bz); (* show sat table *)
    ELSE wgs84(lat, long, heig, bx, by, bz); END;
    satposits(satspos, bx, by, bz, userweek, weekms);
--    dgps(sats, satcnt, satspos, bx, by, bz, mylat, mylong, myhigh, rangcor);
    speedrange(sats, satcnt, satspos);

    speed:=0.0;
    dir:=0.0;
    climb:=0.0;
    speeds:=0;
    FOR i:=0 TO tries-1 DO
      IF stats[i].res*SET16{3..4}=SET16{} THEN
        cnt:=0;
        FOR bit:=0 TO satcnt-1 DO
          IF (bit IN stats[i].satset) & NOT sats[bit].badspeed THEN
            satspos[sats[bit].almidx].range:=LFLOAT(VAL(CARDINAL,MAX(INTEGER)
            - sats[bit].rang))*PRM + sats[bit].userspeed;
            satnums[cnt]:=sats[bit].almidx;
            INC(cnt);
          END;
        END;
        IF (cnt=4) & (get4sats(satspos, satnums, 0, 0.0, vlat, vlong, vheig)>0) THEN
          speed:=speed+neardist(stats[i].lat, stats[i].long, vlat, vlong);
          dirmed(dir, azimuth(vlat, vlong, stats[i].lat, stats[i].long), speeds);
          climb:=climb+stats[i].heig-vheig;
          INC(speeds);
--WrInt(i,2); WrFixed((stats[i].lat-vlat)/RAD*111111.0, 1,6);
--WrFixed((stats[i].long-vlong)/RAD*111111.0, 1,6);
--WrFixed(stats[i].heig-vheig, 1,9);  WrStrLn("");

        END;
      END;
    END;
    IF speeds>0 THEN
      speed:=speed/LFLOAT(speeds);
      dir:=dir/LFLOAT(speeds);
      climb:=climb/LFLOAT(speeds);
    END;
(*
WrFixed(speed, 1,6);
WrFixed(dir, 1,6);
WrStrLn("=speed/dir");
*)
-----speed


    min:=0.0;
    WrStrLn("prn az    ele       clc      dopp   freq0");
    FOR i:=1 TO satcnt DO
      j:=sats[i-1].almidx;
      IF j<=HIGH(satspos) THEN
        IF NOT sats[i-1].badspeed THEN INC(goodsats) END;
        WrInt(sats[i-1].prn, 2);
        WrFixed(satspos[j].azimuth/RAD, 1, 7);
        WrFixed(satspos[j].elevation/RAD, 1, 6); 
        WrFixed(satspos[j].clk, 0, 10);
        WrFixed(satspos[j].doppler, 1, 10);
        WrInt(sats[i-1].rang1, 8);
--        WrFixed(LFLOAT(sats[i-1].freq0)*(WAVLEN/256.0), 3, 10);
        aa:=LFLOAT(sats[i-1].rang1)*(WAVLEN/256.0)-satspos[j].doppler;
        WrFixed(aa, 3, 10);
--      WrFixed(ABS(aa-min), 3, 10);
        min:=aa;
        IF sats[i-1].badspeed THEN WrStr(" v") END;
        WrStrLn("");
      END;
    END;
  END;


  RETURN ret

END getposit;


PROCEDURE readalmanach(fnsem, fnyuma, fnrinex:ARRAY OF CHAR; secondinweek:CARDINAL;
                       VAR tilltime:CARDINAL; verb:BOOLEAN):BOOLEAN;
CONST MAXTRUST=6*3600; (* seconds from oldest to newest entry *)
      MINDATE=(2010-1970)*3600*24*365;
TYPE
    SEM_structAlmanac=RECORD
      toa:CARDINAL;          (* almanac time of applicability (reference time [s]*)
      week,                  (* 10 bit gps week 0-1023 (user must account for week rollover) [week] *)   
      prn,                   (* GPS prn number *)  
      reserved,              (* reserved *)
      svn:CARD16;            (* Satellite vehicle number *)
      ura,                   (* User Range Accuracy lookup code, [0-15], see p. 83 GPSICD200C, 0 is excellent, 15 is use at own risk *)
      health,                (* 0=healthy, unhealthy otherwise  [], subframe 4 and 5, page 25 six-bit health code *)
      config_code:CARD8;      (* configuration code   [], if >=9 Anti-Spoofing is on *)
  
  (* this inicator is not part of the SEM standard but is added by the user if known *)
      is_af0_af1_high_precision:CHAR;  (* indicates precision of af0 and af1 [1=high precision,0=low precision] (22&16 bits, ephemeris source) vs (11&11 bits, almanac source), 0 is typical for most SEM sources *)

      ecc,          (* eccentricity *)
      i0,           (* orbital inclination at reference time                        [rad] *)
      omegadot,     (* rate of right ascension                                      [rad/s] *)
      sqrta,        (* square root of the semi-major axis                           [m^(1/2)]*)
      omega0,       (* longitude of ascending node of orbit plane at weekly epoch   [rad] *)
      w,            (* argument of perigee                                          [rad] *)
      m0,           (* mean anomaly at reference time                               [rad] *)
      af0,          (* polynomial clock correction coefficient (clock bias)         [s],   Note: parameters from ephemeris preferred vs almanac (22 vs 11 bits) *)
      af1:LONGREAL; (* polynomial clock correction coefficient (clock drift)        [s/s], Note: parameters from ephemeris preferred vs almanac (16 vs 11 bits) *)
    END;

    YUMA_structAlmanac=RECORD
      reserved1,
      week,                  (* 10 bit gps week 0-1023 (user must account for week rollover) [week] *)   
      prn:CARD16;            (* GPS prn number *)  
      health:CARD8;          (* 0=healthy, unhealthy otherwise  [], subframe 4 and 5, page 25 six-bit health code *)
      ecc,          (* eccentricity *)
      toa,          (* time of applicability *)
      i0,           (* orbital inclination at reference time                        [rad] *)
      omegadot,     (* rate of right ascension                                      [rad/s] *)
      sqrta,        (* square root of the semi-major axis                           [m^(1/2)]*)
      omega0,       (* longitude of ascending node of orbit plane at weekly epoch   [rad] *)
      w,            (* argument of perigee                                          [rad] *)
      m0,           (* mean anomaly at reference time                               [rad] *)
      af0,          (* polynomial clock correction coefficient (clock bias)         [s],   Note: parameters from ephemeris preferred vs almanac (22 vs 11 bits) *)
      af1:LONGREAL; (* polynomial clock correction coefficient (clock drift)        [s/s], Note: parameters from ephemeris preferred vs almanac (16 vs 11 bits) *)
    END;

VAR fd:INTEGER;
    cnt:CARD8;
    ok:BOOLEAN;
    i, j, ri, ti, minti, nearmed:CARDINAL;
    alm:ARRAY[0..31] OF SEM_structAlmanac;
    yumaalm:ARRAY[0..31] OF YUMA_structAlmanac;
    rinexalm:ARRAY[0..32*96-1] OF structEphemeris;
    min:ARRAY[0..31] OF CARDINAL;
BEGIN
  FILL(ADR(alm), 0C, SIZE(alm));
  FILL(ADR(calm), 0C, SIZE(calm));
  FILL(ADR(yumaalm), 0C, SIZE(yumaalm));
  FILL(ADR(rinexalm), 0C, SIZE(rinexalm));
  semok:=(fnsem[0]<>0C) & (SEM_ReadAlmanacDataFromFile(fnsem, ADR(alm), HIGH(alm)+1, cnt));
  IF semok THEN
    FOR i:=0 TO HIGH(alm) DO alm[i].i0:=alm[i].i0 + 0.3*PI END;
  END;
  yumaok:=(fnyuma[0]<>0C) & YUMA_ReadAlmanacDataFromFile(fnyuma, ADR(yumaalm), HIGH(yumaalm)+1, cnt);
  rinexok:=(fnrinex[0]<>0C) & RINEX_DecodeGPSNavigationFile(fnrinex, ADR(rinexklobuchar),
                              ADR(rinexalm), HIGH(rinexalm), ri);

  IF rinexok & (ri>0) THEN
    IF verb THEN WrInt(ri, 1);WrStrLn("=rec") END;
    FOR i:=0 TO HIGH(min) DO min[i]:=0 END;
    ti:=(time()-7200+3600*24*4) MOD 604800;
    IF verb THEN  WrInt(ti, 12); WrInt(secondinweek, 12); WrStrLn("=ti secondinweek") END;
    FOR j:=0 TO ri-1 DO
      i:=rinexalm[j].prn;
      ti:=(rinexalm[j].tow+SECONDSINWEEK-secondinweek) MOD SECONDSINWEEK;
      IF (i>0) & (i<=HIGH(min)+1) & (ti>SECONDSINWEEK DIV 2) & (min[i-1]<ti) THEN

        DEC(i);
        min[i]:=ti;

--WrInt(i, 8); WrInt(rinexalm[j].week, 8); WrInt(rinexalm[j].tow_week, 8); WrInt(rinexalm[j].tow, 8);
--wrdate(VAL(CARDINAL,rinexalm[j].tow)+VAL(CARDINAL,rinexalm[j].week)*(7*3600*24)+315964800);
--WrStrLn(" ---prn tow  week");

        calm[i].week    :=rinexalm[j].week;
        calm[i].prn     :=rinexalm[j].prn;
        calm[i].health  :=VAL(CARD8, rinexalm[j].health);
        calm[i].ecc     :=rinexalm[j].ecc;
        calm[i].i0      :=rinexalm[j].i0;
        calm[i].omegadot:=rinexalm[j].omegadot;
        calm[i].sqrta   :=rinexalm[j].sqrta;
        calm[i].omega0  :=rinexalm[j].omega0;
        calm[i].w       :=rinexalm[j].w;
        calm[i].m0      :=rinexalm[j].m0;
        calm[i].af0     :=rinexalm[j].af0;
        calm[i].af1     :=rinexalm[j].af1;
        calm[i].af2     :=rinexalm[j].af2;
        calm[i].toe     :=rinexalm[j].toe;
        calm[i].toc     :=rinexalm[j].toc;
        calm[i].tow_week:=rinexalm[j].tow_week;
        calm[i].tow     :=rinexalm[j].tow;
        calm[i].tgd     :=rinexalm[j].tgd;
        calm[i].delta_n :=rinexalm[j].delta_n;
        calm[i].idot    :=rinexalm[j].idot;
        calm[i].cuc     :=rinexalm[j].cuc;
        calm[i].cus     :=rinexalm[j].cus;
        calm[i].crc     :=rinexalm[j].crc;
        calm[i].crs     :=rinexalm[j].crs;
        calm[i].cic     :=rinexalm[j].cic;
        calm[i].cis     :=rinexalm[j].cis;

        calm[i].treal:=VAL(CARDINAL,calm[i].tow)+VAL(CARDINAL,calm[i].week)*(7*3600*24)+315964800;
      END;
    END;

    j:=0;
    nearmed:=MAX(INTEGER);
    REPEAT                                            (* find median time *)
      ti:=0;      
      FOR i:=0 TO HIGH(calm) DO
        IF calm[j].treal<calm[i].treal THEN INC(ti) END;
      END;
      ti:=ABS(VAL(INTEGER,HIGH(calm) DIV 2)-VAL(INTEGER,ti)); 
      IF ti<nearmed THEN nearmed:=ti; tilltime:=calm[j].treal END;
--IF verb THEN WrStr("sat time:"); WrInt(j, 8); WrInt(ti, 8); WrInt(nearmed, 8); WrStr(" "); wrdate(calm[j].treal); WrStrLn("") END;
      INC(j);
    UNTIL j>HIGH(calm); 

    IF verb THEN WrStr("median last almanach time:"); wrdate(tilltime); WrStr(" ") END;
 

(*
    minti:=MAX(CARDINAL);                     (* find oldest entry for difference to newest *)
    FOR i:=0 TO 31 DO
      ti:=VAL(CARDINAL,calm[i].tow)+VAL(CARDINAL,calm[i].week)*(7*3600*24)+315964800;
      IF (ti>MINDATE) & (ti<minti) THEN minti:=ti END;              (* oldest entry *)
    END;

    minti2:=MAX(CARDINAL);                     (* find second oldest entry for difference to newest *)
    FOR i:=0 TO 31 DO
      ti:=VAL(CARDINAL,calm[i].tow)+VAL(CARDINAL,calm[i].week)*(7*3600*24)+315964800;
      IF (ti>minti) & (ti<minti2) THEN minti2:=ti END;              (* second oldest entry *)
    END;
    minti:=minti2;

    tilltime:=0;
    FOR i:=0 TO HIGH(calm) DO
-- WrInt(calm[i].tow, 12); WrInt(calm[i].week, 10);
      ti:=VAL(CARDINAL,calm[i].tow)+VAL(CARDINAL,calm[i].week)*(7*3600*24)+315964800;   
      IF (ti<minti+MAXTRUST) & (ti>tilltime) THEN tilltime:=ti END;  (* newest trusted entry as hint for alm timeout *)
      IF verb THEN 
        wrdate(ti);
        IF i MOD 4=3 THEN WrStrLn("") ELSE WrStr(" ") END;
      END;
    END;
*)

  ELSIF semok THEN
    FOR i:=0 TO HIGH(alm) DO
      calm[i].toa     :=alm[i].toa;
      calm[i].week    :=alm[i].week;
      calm[i].prn     :=alm[i].prn;
      calm[i].health  :=alm[i].health;
      calm[i].ecc     :=alm[i].ecc;
      calm[i].i0      :=alm[i].i0;
      calm[i].omegadot:=alm[i].omegadot;
      calm[i].sqrta   :=alm[i].sqrta;
      calm[i].omega0  :=alm[i].omega0;
      calm[i].w       :=alm[i].w;
      calm[i].m0      :=alm[i].m0;
      calm[i].af0     :=alm[i].af0;
      calm[i].af1     :=alm[i].af1;
    END;


  ELSIF yumaok THEN
--    IF verb THEN
(*
      FOR i:=0 TO HIGH(calm) DO
        j:=0;
        WHILE (j<=HIGH(yumaalm)) & (alm[i].prn<>yumaalm[j].prn) DO INC(j) END;
        IF j<=HIGH(yumaalm) THEN
          WrInt(alm[i].prn, 3); WrStrLn("");
          WrInt(alm[i].toa, 25); WrInt(TRUNC(yumaalm[j].toa), 25); WrStrLn("");
          WrInt(alm[i].week, 25); WrInt(yumaalm[j].week, 25); WrStrLn("");
          WrInt(alm[i].health, 25); WrInt(yumaalm[j].health, 25); WrStrLn("");
          WrFixed(alm[i].ecc, 12, 25); WrFixed(yumaalm[j].ecc, 12, 25); WrStrLn("");
          WrFixed(alm[i].i0, 12, 25);  WrFixed(yumaalm[j].i0, 12, 25); WrStrLn("");
          WrFixed(alm[i].omegadot, 12, 25); WrFixed(yumaalm[j].omegadot, 12, 25); WrStrLn("");
          WrFixed(alm[i].sqrta, 12, 25); WrFixed(yumaalm[j].sqrta, 12, 25); WrStrLn("");
          WrFixed(alm[i].omega0, 12, 25); WrFixed(yumaalm[j].omega0, 12, 25); WrStrLn("");
          WrFixed(alm[i].w, 12, 25); WrFixed(yumaalm[j].w, 12, 25); WrStrLn("");
          WrFixed(alm[i].m0, 12, 25); WrFixed(yumaalm[j].m0, 12, 25); WrStrLn("");
          WrFixed(alm[i].af0, 12, 25); WrFixed(yumaalm[j].af0, 12, 25); WrStrLn("");
          WrFixed(alm[i].af1, 12, 25); WrFixed(yumaalm[j].af1, 12, 25); WrStrLn("");
        END;
*)
--      END;

    FOR i:=0 TO HIGH(calm) DO
      calm[i].toa     :=TRUNC(yumaalm[i].toa);
      calm[i].week    :=yumaalm[i].week;
      calm[i].prn     :=yumaalm[i].prn;
      calm[i].health  :=yumaalm[i].health;
      calm[i].ecc     :=yumaalm[i].ecc;
      calm[i].i0      :=yumaalm[i].i0;
      calm[i].omegadot:=yumaalm[i].omegadot;
      calm[i].sqrta   :=yumaalm[i].sqrta;
      calm[i].omega0  :=yumaalm[i].omega0;
      calm[i].w       :=yumaalm[i].w;
      calm[i].m0      :=yumaalm[i].m0;
      calm[i].af0     :=yumaalm[i].af0;
      calm[i].af1     :=yumaalm[i].af1;
    END;
  END;

  RETURN semok OR yumaok OR rinexok
END readalmanach;


END gpspos.
(*
    toa:CARDINAL;          (* almanac time of applicability (reference time [s]*)
    week,                  (* 10 bit gps week 0-1023 (user must account for week rollover) [week] *
    prn,                   (* GPS prn number *)
    reserved,              (* reserved *)
    svn:CARD16;            (* Satellite vehicle number *)
    ura,                   (* User Range Accuracy lookup code, [0-15], see p. 83 GPSICD200C, 0 is e
    health,                (* 0=healthy, unhealthy otherwise  [], subframe 4 and 5, page 25 six-bit
    config_code:CARD8;      (* configuration code   [], if >=9 Anti-Spoofing is on *)

  (* this inicator is not part of the SEM standard but is added by the user if known *)
    is_af0_af1_high_precision:CHAR;  (* indicates precision of af0 and af1 [1=high precision,0=low

    ecc,          (* eccentricity *)
    i0,           (* orbital inclination at reference time                        [rad] *)
    omegadot,     (* rate of right ascension                                      [rad/s] *)
    sqrta,        (* square root of the semi-major axis                           [m^(1/2)]*)
    omega0,       (* longitude of ascending node of orbit plane at weekly epoch   [rad] *)
    w,            (* argument of perigee                                          [rad] *)
    m0,           (* mean anomaly at reference time                               [rad] *)
    af0,          (* polynomial clock correction coefficient (clock bias)         [s],   Note: para
    af1:LONGREAL; (* polynomial clock correction coefficient (clock drift)        [s/s], Note: para


    reserved1,
    week,                  (* 10 bit gps week 0-1023 (user must account for week rollover) [week] *
    prn:CARD16;            (* GPS prn number *)
    health:CARD8;          (* 0=healthy, unhealthy otherwise  [], subframe 4 and 5, page 25 six-bit
    ecc,          (* eccentricity *)
    toa,          (* time of applicability *)
    i0,           (* orbital inclination at reference time                        [rad] *)
    omegadot,     (* rate of right ascension                                      [rad/s] *)
    sqrta,        (* square root of the semi-major axis                           [m^(1/2)]*)
    omega0,       (* longitude of ascending node of orbit plane at weekly epoch   [rad] *)
    w,            (* argument of perigee                                          [rad] *)
    m0,           (* mean anomaly at reference time                               [rad] *)
    af0,          (* polynomial clock correction coefficient (clock bias)         [s],   Note: para
    af1:LONGREAL; (* polynomial clock correction coefficient (clock drift)        [s/s], Note: para


6246 795 1304 doppler lock error

0kmh (not im mean)
max dop delete dynamic
lim age of almanach
--vrms limit
if wrong alt wait a time then send without alt
*)