fork download
copy 
  1. program stdalone;
  2.  
  3. type
  4.    mat96 = array[1..32,1..3] of real;
  5.    vecb32 = array[1..32] of boolean;
  6.    vec32 = array[1..32] of real;
  7.    vec16 = array[1..16] of real;
  8.    vec8 = array[1..8] of real;
  9.    vec3 = array[1..3] of real;
  10.    mat16 = array[1..4,1..4] of real;
  11.  
  12. const
  13.    We = 7.292115E-5;            {WGS-84 earth rotation rate}
  14.    c = 299792458.0;             {WGS-84 speed of light}
  15.    pi = 3.1415926535898         {WGS 84 value of pi};
  16.    Wedot = 7.2921151467E-5;     {WGS 84 value of earth's rotation rate}
  17.    mu = 3.986005E+14;           {WGS 84 value of earth's univ. grav. par.}
  18.    F = -4.442807633E-10;        {relativistic correction term constant}
  19.    a = 6378137.0;               {WGS-84 earth's semi major axis}
  20.    b = 6356752.31;              {WGS-84 earth's semi minor axis}
  21.    e1sqr = (a*a - b*b) / (a*a); {first  numerical eccentricity}
  22.    e2sqr = (a*a - b*b) / (b*b); {second numerical eccentricity}
  23.  
  24. var
  25.    inp, out: text;
  26.    Ttr, tau, Trc, T, Trel: real;
  27.    Az, El, Cr, alpha, dTclck, dTiono, dRtrop, Lat, Lon: real;
  28.    prn, i: integer;
  29.    eph: array[1..32,1..16] of real;
  30.    clk: array[1..32,1..5] of real;
  31.    ion: vec8;
  32.    Praw, Pcor: vec32;
  33.    Xs: mat96;
  34.    SV: vecb32;
  35.    P: vec32;
  36.    Xr, Xlla, tmp3: vec3;
  37.    status: boolean;
  38.    tmp16: vec16;
  39.    tmp5: array[1..5] of real;
  40.  
  41. {***************************************************************************}
  42.  
  43. procedure LLA2XYZ(Xi : vec3;      {lat [rad], lon [rad] alt [m]}
  44.               var Xo : vec3);     {ECEF X [m], Y [m] , Z[m]}
  45. {this procedure converts WGS-84 Lat, Lon and Alt [above ellipsoid]
  46.  to ECEF XYZ}
  47.  
  48. var
  49.    N: real;
  50. begin
  51. N := a / sqrt(1.0 - e1sqr * sin(Xi[1]) * sin(Xi[1]));
  52. Xo[1] := (N                 + Xi[3]) * cos(Xi[1]) * cos(Xi[2]);
  53. Xo[2] := (N                 + Xi[3]) * cos(Xi[1]) * sin(Xi[2]);
  54. Xo[3] := (N * (1.0 - e1sqr) + Xi[3]) * sin(Xi[1]);
  55. end; {procedure LLA2XYZ}
  56.  
  57. {***************************************************************************}
  58.  
  59. procedure XYZ2LLA(Xi : vec3;      {ECEF X [m], Y [m] , Z[m]}
  60.               var Xo : vec3);     {lat [rad], lon [rad] alt [m]}
  61. {this procedure converts WGS-84 ECEF XYZ to Lat, Lon, Alt [above ellipsoid]}
  62.  
  63. var
  64.    p, T, sT, cT, N, sig: real;
  65. begin
  66. p := sqrt(Xi[1] * Xi[1] + Xi[2] * Xi[2]);
  67. T := arctan((Xi[3] * a) / (p * b));
  68. sT := sin(T); cT := cos(T);
  69. Xo[1] := arctan((Xi[3] + e2sqr * b * sT * sT * sT)
  70.          / (p - e1sqr * a * cT * cT * cT));
  71. if Xi[2] <> 0.0 then sig := Xi[2] / abs(Xi[2]) else sig := 1.0;
  72. if Xi[1] = 0.0 then Xo[2] := sig * pi / 2.0 else
  73.    begin
  74.    Xo[2] := arctan(Xi[2]/Xi[1]);
  75.    if (Xi[1] < 0.0) and (Xi[2] >= 0.0) then Xo[2] := Xo[2] + pi;
  76.    if (Xi[1] < 0.0) and (Xi[2] < 0.0) then Xo[2] := Xo[2] - pi
  77.    end;
  78. N := a / sqrt(1.0 - e1sqr * sin(Xo[1]) * sin(Xo[1]));
  79. Xo[3] := p / cos(Xo[1]) - N;
  80.  
  81. end; {procedure XYZ2LLA}
  82.  
  83. {***************************************************************************}
  84.  
  85. procedure satpos(eph : vec16;    {ephemeris}
  86.                  Ttr : real;     {satellite GPS time}
  87.              var Trel: real;     {relativistic correction term}
  88.              var X   : vec3);    {satellite position}
  89.  
  90. var
  91.   M0, dn, ec, A, W0, i0, w, Wdot, Cuc, Cus, Crc, Crs, Cic, Cis, Toe, Idot: real;
  92.   T, n0, n, M, E, Eold, snu, cnu, nu, phi, du, dr, di, u, r, i, Xdash, Ydash, Wc: real;
  93.   k: integer;
  94.  
  95. begin
  96. {for clarity of code, copy the ephemeris parameters and convert to radians}
  97. Crs := eph[1];
  98. dn  := eph[2] * pi;
  99. M0  := eph[3] * pi;
  100. Cuc := eph[4];
  101. ec  := eph[5];
  102. Cus := eph[6];
  103. A   := eph[7] * eph[7];
  104. Toe := eph[8];
  105. Cic := eph[9];
  106. W0  := eph[10] * pi;
  107. Cis := eph[11];
  108. i0  := eph[12] * pi;
  109. Crc := eph[13];
  110. w   := eph[14] * pi;
  111. Wdot:= eph[15] * pi;
  112. idot:= eph[16] * pi;
  113.  
  114. T:= Ttr - Toe;
  115. if T >  302400 then T := T - 604800;
  116. if T < -302400 then T := T + 604800;
  117.  
  118. n0 := sqrt(mu / (A*A*A));
  119. n := n0 + dn;
  120.  
  121. M := M0 + n*T;
  122. E := M; {start value for E}
  123. repeat
  124.    Eold := E;
  125.    E := M + ec * sin(E);
  126.    until abs(E - Eold) < 1.0e-8;
  127.  
  128. snu := sqrt(1 - ec*ec) * sin(E) / (1 - ec*cos(E)); 
  129. cnu := (cos(E) - ec) / (1 - ec*cos(E)); 
  130. if cnu = 0 then nu := pi/2 * snu / abs(snu)
  131. else if (snu = 0) and (cnu > 0) then nu := 0
  132. else if (snu = 0) and (cnu < 0) then nu := pi
  133. else nu := arctan(snu/cnu)
  134. + ord(cnu<0) * pi * snu / abs(snu);
  135.  
  136. phi := nu + w;
  137.  
  138. du := Cuc*cos(2*phi) + Cus*sin(2*phi);
  139. dr := Crc*cos(2*phi) + Crs*sin(2*phi);
  140. di := Cic*cos(2*phi) + Cis*sin(2*phi);
  141.  
  142. u := phi + du;
  143. r := A*(1 - ec*cos(E)) + dr;
  144. i := i0 + idot*T +di;
  145.  
  146. Xdash := r*cos(u);
  147. Ydash := r*sin(u);
  148.  
  149. Wc:= W0 + (Wdot - Wedot)*T - Wedot*Toe;
  150.  
  151. X[1] := Xdash*cos(Wc) - Ydash*cos(i)*sin(Wc);
  152. X[2] := Xdash*sin(Wc) + Ydash*cos(i)*cos(Wc);
  153. X[3] := Ydash*sin(i);
  154.  
  155. Trel := F * ec * eph[7] * sin(E) {relativistic correction term}
  156.  
  157. end; {procedure satpos}
  158.  
  159. {***************************************************************************}
  160.  
  161. procedure calcAzEl(Xs,              {satellite ECEF XYZ}
  162.                    Xu  :  vec3;     {user ECEF XYZ}
  163.                var Az,              {azimuth [rad]}
  164.                    El  : real;      {elevation [rad]}
  165.                var stat: boolean);  {calculation succeeded: stat = true}
  166.  
  167. var
  168.    R, p, x, y, z, s: real;
  169.    e: array[1..3,1..3] of real;
  170.    i, k: integer;
  171.    d: vec3;
  172.  
  173. begin
  174.  
  175. x := Xu[1];
  176. y := Xu[2];
  177. z := Xu[3];
  178. p := sqrt(x*x + y*y);
  179. if p = 0.0 then stat := false else stat := true;
  180.  
  181. if stat then
  182.    begin
  183.  
  184.    R := sqrt(x*x + y*y + z*z);
  185.    e[1,1] := - y / p;
  186.    e[1,2] := + x / p;
  187.    e[1,3] := 0.0;
  188.    e[2,1] := - x*z / (p*R);
  189.    e[2,2] := - y*z / (p*R);
  190.    e[2,3] := p / R;
  191.    e[3,1] := x / R;
  192.    e[3,2] := y / R;
  193.    e[3,3] := z / R;
  194.  
  195.    for k := 1 to 3 do
  196.       begin
  197.       d[k] := 0.0;
  198.       for i := 1 to 3 do d[k] := d[k] + (Xs[i] - Xu[i]) * e[k,i]
  199.       end;
  200.  
  201.    s := d[3] / sqrt(d[1]*d[1] + d[2]*d[2] + d[3]*d[3]);
  202.    if s = 1.0 then El := 0.5 * pi else El := arctan(s / sqrt(1.0 - s*s));
  203.  
  204.    if (d[2] = 0.0) and (d[1] > 0.0) then Az := 0.5 * pi else
  205.       if (d[2] = 0.0) and (d[1] < 0.0) then Az := 1.5 * pi else
  206.          begin
  207.          Az := arctan(d[1] / d[2]);
  208.          if d[2] < 0.0 then Az := Az + pi else
  209.             if (d[2] > 0.0) and (d[1] < 0.0) then Az := Az + 2.0 * pi
  210.          end;
  211.  
  212.    end;
  213.  
  214. end; {procedure calcAzEl}
  215.  
  216. {***************************************************************************}
  217.  
  218. procedure ionocorr (ion   :vec8;    {iono correction coefficients from
  219.                                      nav message}
  220.                     Latu,           {user's latitude [rad]}
  221.                     Lonu,           {user's longitude [rad]}
  222.                     Az,             {SV azimuth [rad]}
  223.                     El,             {SV elevation [rad]}
  224.                     Ttr   : real;   {SV signal transmission time [sec]}
  225.                 var dTiono: real);  {Ionospheric delay [sec]}
  226.  
  227. var
  228.   phi, Lati, Loni, Latm, T, F, x, per, amp: real;
  229.   a0, a1, a2, a3, b0, b1, b2, b3: real;
  230.  
  231. begin
  232.  
  233. {for clarity copy array}
  234. a0 := ion[1];
  235. a1 := ion[2];
  236. a2 := ion[3];
  237. a3 := ion[4];
  238. b0 := ion[5];
  239. b1 := ion[6];
  240. b2 := ion[7];
  241. b3 := ion[8];
  242.  
  243. {convert from radians to semicircles}
  244. Latu := Latu / pi; Lonu := Lonu / pi; Az := Az / pi; El := El / pi;
  245.  
  246. {calculation}
  247. phi := 0.0137 / (El + 0.11) - 0.022;
  248. Lati := Latu + phi * cos (Az * pi);
  249. if Lati > +0.416 then Lati := +0.416 else if Lati < -0.416 then Lati := -0.416;
  250. Loni := Lonu + phi * sin(Az * pi) / cos(Lati * pi);
  251. Latm := Lati + 0.064 * cos((Loni - 1.617) * pi);
  252. T := 4.32E+4 * Loni + Ttr;
  253. while T >= 86400 do T := T - 86400 else while T < 0 do T := T + 86400;
  254. F := 1.0 + 16.0 * (0.53 - El) * (0.53 - El) * (0.53 - El);
  255. per := b0 + b1 * Latm + B2 * Latm * Latm + b3 * Latm * Latm * Latm;
  256. if per < 72000.0 then per := 72000.0;
  257. x := 2 * pi * (T - 50400.0) / per;
  258. amp := a0 + a1 * Latm + A2 * Latm * Latm + a3 * Latm * Latm * Latm;
  259. if amp < 0.0 then amp := 0.0;
  260. if abs(x) >= 1.57 then dTiono := F * 5.0E-9 else
  261.    dTiono := F * (5.0E-9 + amp * (1.0 - x * x / 2.0 + x * x * x * x /24.0))
  262.  
  263. end; {procedure ionocorr}
  264.  
  265. {***************************************************************************}
  266.  
  267. {At many places in the following procedure solve the subdeterminant value of a
  268.  4 x 4 array is required. For convenience the function sub is defined below}
  269.  
  270. function sub (A : mat16;      {input 4 x 4 array}
  271.               r,              {row number to be deleted}
  272.               c : integer     {column number to be deleted}
  273.               ) : double;     {value of 3 x 3 subdeterminant}
  274.  
  275. var
  276.   B: array[1..3,1..3] of double;
  277.   i, i1, j, j1: integer;
  278.  
  279. begin
  280. i1 := 0;
  281. for i := 1 to 4 do if i <> r then
  282.    begin
  283.    i1 := i1 + 1;
  284.    j1 := 0;
  285.    for j := 1 to 4 do if j <> c then
  286.       begin
  287.       j1 := j1 + 1;
  288.       B[i1,j1] := A[i,j]
  289.       end
  290.    end;
  291. sub := + B[1,1]*(B[2,2]*B[3,3] - B[2,3]*B[3,2])
  292.        - B[2,1]*(B[1,2]*B[3,3] - B[3,2]*B[1,3])
  293.        + B[3,1]*(B[1,2]*B[2,3] - B[1,3]*B[2,2]);
  294. end; {function sub}
  295.  
  296. {***************************************************************************}
  297.  
  298. procedure solve(Xs    : mat96;    {array with 3 columns and 32 rows
  299.                                    for the coordinates of the sat's}
  300.                 SV    : vecb32;   {valid prn's}
  301.                 P     : vec32;    {pseudoranges}
  302.             var Xr    : vec3;     {input of initial guess, output of
  303.                                    final position}
  304.             var Cr    : real;     {receiver clock error}
  305.             var status: boolean); {true: calculation OK, false: no solution}
  306.  
  307. {procedure solve requires the following types to be declared in the
  308.  main body of the program:
  309.  type
  310.     mat96 = array[1..32,1..3] of real;
  311.     vecb32 = array[1..32] of boolean;
  312.     vec32 = array[1..32] of real;
  313.     vec3 = array[1..3] of real;
  314.     mat16 = array[1..4,1..4] of real;}
  315.  
  316. var
  317.    prn, it, i, j, k: integer;
  318.    R, L: array[1..32] of real;
  319.    A: array[1..32,1..4] of real;
  320.    AL: array[1..4] of real;
  321.    AA, AAi: mat16;
  322.    n: longint;
  323.    det: real;
  324.    D: array[1..4] of real;
  325.  
  326.  
  327. begin {procedure solve}
  328.  
  329. it := 0; {iteration counter}
  330.  
  331. repeat {iterations}
  332.  
  333.    it :=it + 1; {increase iteration counter}
  334.  
  335.    for prn := 1 to 32 do if SV[prn] then
  336.       begin
  337.       R[prn] :=  {range from receiver to satellite}
  338.          sqrt((Xr[1] - Xs[prn,1]) * (Xr[1] - Xs[prn,1])
  339.             + (Xr[2] - Xs[prn,2]) * (Xr[2] - Xs[prn,2])
  340.             + (Xr[3] - Xs[prn,3]) * (Xr[3] - Xs[prn,3]));
  341.       L[prn] := P[prn] - R[prn]; {range residual value}
  342.       for k := 1 to 3 do A[prn,k] := (Xr[k] - Xs[prn,k]) / R[prn];
  343.       A[prn,4] := -1.0 {A is the geometry matrix or model matrix}
  344.       end;
  345.  
  346.    For k :=1 to 4 do {calculate A.L}
  347.       begin
  348.       AL[k] := 0.0;
  349.       for prn := 1 to 32 do if SV[prn] then
  350.          AL[k] := AL[k] + A[prn,k] * L[prn]
  351.       end;
  352.  
  353.    for k := 1 to 4 do for i := 1 to 4 do {calculate A.A}
  354.       begin
  355.       AA[k,i] :=0.0;
  356.       for prn := 1 to 32 do if SV[prn] then
  357.          AA[k,i] := AA[k,i] + A[prn,k] * A[prn,i]
  358.       end;
  359.  
  360.    {invert A.A}
  361.    det := + AA[1,1] * sub(AA,1,1) - AA[2,1] * sub(AA,2,1)
  362.           + AA[3,1] * sub(AA,3,1) - AA[4,1] * sub(AA,4,1);
  363.    if det = 0.0 then status := false else
  364.       begin
  365.  
  366.       status := true;
  367.  
  368.       for k := 1 to 4 do for i := 1 to 4 do
  369.          begin
  370.          n:= k + i; if odd(n) then j := -1 else j :=1;
  371.          AAi[k,i] := j * sub(AA,i,k) / det
  372.          end;
  373.  
  374.       {calculate (invA.A).(A.L)}
  375.       for k := 1 to 4 do
  376.          begin
  377.          D[k] := 0.0;
  378.          for i := 1 to 4 do D[k] := D[k] + AAi[k,i] * AL[i]
  379.          end;
  380.  
  381.       {update position}
  382.       for k := 1 to 3 do Xr[k] := Xr[k] + D[k];
  383.  
  384.       end;
  385.  
  386.    until (it = 6){there is something wrong if more than 6 iterations are required}
  387.       or ((abs(D[1]) + abs(D[2]) + abs(D[3])) < 1.0E-2) {iteration criterion}
  388.          or (not(status)); {calculation not succeeded}
  389.  
  390. Cr := D[4]; {receiver clock error}
  391.  
  392. if it = 6 then begin writeln('solve it : ',it); status := false end; {iteration not succeeded}
  393.  
  394. end; {procedure solve}
  395.  
  396. {***************************************************************************}
  397.  
  398. begin {main}
  399. {the following data should be available:
  400.  1. Pseudorange with receiver time of reception for each SV
  401.  2. Ephemeris and almanac for each SV
  402.  3. Iono coefficients}
  403.  
  404. {open input datafile}
  405. assign(inp,'inp.txt'); reset(inp);
  406.  
  407. {read GPS time of reception}
  408. readln(inp); {skip comment line}
  409. readln(inp,Trc);
  410.  
  411. {read iono coefficients}
  412. readln(inp); {skip comment line}
  413. for i := 1 to 8 do readln(inp,ion[i]);
  414.  
  415. readln(inp); {skip comment line}
  416. {read pseudoranges}
  417. for prn := 1 to 32 do SV[prn] := false;
  418. repeat
  419.    read(inp,prn);
  420.    if prn <> 0 then
  421.       begin
  422.       readln(inp,Praw[prn]);
  423.       SV[prn] := true
  424.       end
  425.    else readln(inp);
  426.    until prn = 0;
  427.  
  428. readln(inp); {skip comment line}
  429. {read ephemeris- and clock data}
  430. repeat
  431.    readln(inp,prn);
  432.    for i := 1 to 16 do readln(inp,eph[prn,i]);
  433.    for i := 1 to 5 do readln(inp,clk[prn,i]);
  434.    until eof(inp);
  435. close(inp);
  436.  
  437. {user input of start position}
  438. write('Start position Lat [deg.dec], Lon [deg.dec], Alt [m] : ');
  439.    readln(Xlla[1],Xlla[2],Xlla[3]);
  440. Xlla[1] := Xlla[1] * pi / 180.0; Xlla[2] :=Xlla[2] * pi / 180.0;
  441. {convert lat, ln, alt to ECEF X, Y, Z}
  442. LLA2XYZ(Xlla,Xr);
  443.  
  444. {open output data file}
  445. assign(out,'output.txt'); rewrite(out);
  446.  
  447. {assuming the receiver clock error and initial position not sufficiently
  448.  good known, I make two passes through the processing steps}
  449.  
  450. {PASS 1}
  451. writeln(out,'PASS 1');
  452.  
  453. for prn := 1 to 32 do if SV[prn] then begin {do for each SV}
  454.  
  455.    {set all transit times to nominal value and calculate time of transmission}
  456.    tau := 0.075;
  457.    Ttr := Trc - tau;
  458.  
  459.    {calculate SV position and correct for earth rotation}
  460.    for i := 1 to 16 do tmp16[i] := eph[prn,i];
  461.    satpos(tmp16,Ttr,Trel,tmp3);
  462.    alpha := tau * We;
  463.    Xs[prn,1] := + tmp3[1] * cos(alpha) + tmp3[2] * sin(alpha);
  464.    Xs[prn,2] := - tmp3[1] * sin(alpha) + tmp3[2] * cos(alpha);
  465.    Xs[prn,3] := + tmp3[3];
  466.    writeln(out,'SV     : ',prn:2,Xs[prn,1]:15:3,Xs[prn,2]:15:3,Xs[prn,3]:15:3);
  467.  
  468.    {calculate azimuth and elevation}
  469.    for i := 1 to 3 do tmp3[i] := Xs[prn,i];
  470.    calcAzEl(tmp3,Xr,Az,El,status);
  471.    if not status then
  472.       begin writeln('Error in calcAzEl - check input data'); exit end;
  473.    writeln(out,'Az, El : ',prn:2,Az*180.0/pi:11:3,El*180.0/pi:10:3);
  474.  
  475.    {calculate pseudorange corrections and apply to pseudoranges}
  476.  
  477.    {clock correction}
  478.    T := Ttr - clk[prn,2];
  479.    {correct for week crossover}
  480.    if T >  302400 then T := T - 604800;
  481.    if T < -302400 then T := T + 604800;
  482.  
  483.    dTclck := + clk[prn,5] + clk[prn,4] * T + clk[prn,3] * T * T
  484.              + Trel - clk[prn,1];
  485.  
  486.    {iono correction}
  487.    Lat := Xlla[1] ; Lon := Xlla[2];
  488.    ionocorr(ion,Lat,Lon,Az,El,Ttr,dTiono);
  489.  
  490.    {tropo correction using standard atmosphere values}
  491.    dRtrop := + 2.312 / sin(sqrt(El * El + 1.904E-3))
  492.              + 0.084 / sin(sqrt(El * El + 0.6854E-3));
  493.  
  494.    writeln(out,'Corr   : ',prn:2,dTclck*c:11:3,dTiono*c:10:3,dRtrop:10:3);
  495.  
  496.    {correct pseudorange}
  497.    Pcor[prn] := Praw[prn] + dTclck * c - dTiono * c - dRtrop
  498.  
  499.    end; {do for each SV}
  500.  
  501. {calculate receiver position}
  502. solve(Xs,SV,Pcor,Xr,Cr,status);
  503. if not status then
  504.    begin writeln('Error in solve - check input data'); exit end;
  505. writeln(out,'Pos XYZ: ',Xr[1]:12:3,Xr[2]:12:3,Xr[3]:12:3,Cr:12:3);
  506. {convert back to Lat, Lon, Alt}
  507. XYZ2LLA(Xr,Xlla);
  508.  
  509. {PASS 2 - The receiver position and -clock error is now well enough known
  510.  to calculate the final pseudorange corrections}
  511. writeln(out); writeln(out,'PASS 2');
  512.  
  513. {correct receiver clock}
  514. Trc := Trc + Cr / c;
  515.  
  516. for prn := 1 to 32 do if SV[prn] then begin {do for each SV}
  517.  
  518.    {recalculate transit time and time of transmission}
  519.    tau := (Pcor[prn] + Cr) / c;
  520.    Ttr := Trc - tau;
  521.  
  522.    {recalculate SV position and correct for earth rotation}
  523.    for i := 1 to 16 do tmp16[i] := eph[prn,i];
  524.    satpos(tmp16,Ttr,Trel,tmp3);
  525.    alpha := tau * We;
  526.    Xs[prn,1] := + tmp3[1] * cos(alpha) + tmp3[2] * sin(alpha);
  527.    Xs[prn,2] := - tmp3[1] * sin(alpha) + tmp3[2] * cos(alpha);
  528.    Xs[prn,3] := + tmp3[3];
  529.    writeln(out,'SV     : ',prn:2,Xs[prn,1]:15:3,Xs[prn,2]:15:3,Xs[prn,3]:15:3);
  530.  
  531.    {recalculate azimuth and elevation}
  532.    for i := 1 to 3 do tmp3[i] := Xs[prn,i];
  533.    calcAzEl(tmp3,Xr,Az,El,status);
  534.    if not status then
  535.       begin writeln('Error in calcAzEl - check input data'); exit end;
  536.    writeln(out,'Az, El : ',prn:2,Az*180.0/pi:11:3,El*180.0/pi:10:3);
  537.  
  538.    {recalculate pseudorange corrections and apply to pseudoranges}
  539.  
  540.    {clock correction}
  541.    T := Ttr - clk[prn,2];
  542.    {correct for week crossover}
  543.    if T >  302400 then T := T - 604800;
  544.    if T < -302400 then T := T + 604800;
  545.    dTclck := + clk[prn,5] + clk[prn,4] * T + clk[prn,3] * T * T
  546.              + Trel - clk[prn,1];
  547.  
  548.    {iono correction}
  549.    Lat := Xlla[1]; Lon := Xlla[2];
  550.    ionocorr(ion,Lat,Lon,Az,El,Ttr,dTiono);
  551.  
  552.    {tropo correction using standard atmosphere values}
  553.    dRtrop := + 2.312 / sin(sqrt(El * El + 1.904E-3))
  554.              + 0.084 / sin(sqrt(El * El + 0.6854E-3));
  555.  
  556.    writeln(out,'Corr   : ',prn:2,dTclck*c:11:3,dTiono*c:10:3,dRtrop:10:3);
  557.  
  558.    {correct pseudorange}
  559.    Pcor[prn] := Praw[prn] + dTclck * c - dTiono * c - dRtrop + Cr
  560.  
  561.    end; {do for each SV}
  562.  
  563. {calculate receiver position}
  564. solve(Xs,SV,Pcor,Xr,Cr,status);
  565. if not status then
  566.    begin writeln('Error in solve - check input data'); exit end;
  567. writeln(out,'Pos XYZ: ',Xr[1]:12:3,Xr[2]:12:3,Xr[3]:12:3,Cr:12:3);
  568. {convert back to Lat, Lon, Alt}
  569. XYZ2LLA(Xr,Xlla);
  570. writeln(out,'Pos LLA: ',Xlla[1]*180.0/pi:15:8,Xlla[2]*180.0/pi:15:8,Xlla[3]:12:3);
  571.  
  572. close(out)
  573. end. {main}
  574.  
  575.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
compilation info 
Free Pascal Compiler version 2.2.0 [2009/11/16] for i386
Copyright (c) 1993-2007 by Florian Klaempfl
Target OS: Linux for i386
Compiling prog.pas
prog.pas(93,3) Note: Local variable "k" not used
prog.pas(253,36) Fatal: Syntax error, ";" expected but "ELSE" found
Fatal: Compilation aborted
Error: /usr/bin/ppc386 returned an error exitcode (normal if you did not specify a source file to be compiled)
stdout
copy 
Standard output is empty
https://ideone.com/gJUmb
language:
Pascal (fpc 3.0.4)
created:
14 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!