home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / worldmap / mapvu20.arc / MAPPROJ.PAS < prev    next >
Pascal/Delphi Source File  |  1989-01-11  |  33KB  |  829 lines
  1. Unit mapproj;  { projection stuff for MapView }
  2.  
  3. { Copyright A.J. van den Bogert and Gisbert W.Selke Dec 1988                 }
  4.  
  5. {$UNDEF  DEBUG }                       { DEFINE while debugging }
  6.  
  7. {$IFDEF DEBUG }
  8. {$A+,R+,S+,I+,D+,F-,V-,B-,L+ }      { turn checking on while debugging }
  9. {$ELSE }
  10. {$A+,R-,S-,I+,D-,F-,V-,B-,L+ }
  11. {$ENDIF }
  12. {$M 65500,65500,560000}
  13.  
  14. {$IFDEF CPU87 }
  15. {$N+ }
  16. {$ELSE }
  17. {$N- }
  18. {$ENDIF }
  19.  
  20. Interface
  21.  
  22. Uses  Crt, Graph, mapgraph;
  23.  
  24. Const maxbuf        = 100;      { max. size of curve buffer }
  25.       trigtablesize = 9000;     { size of table for precomputed trig values }
  26.       epsilon       = 1e-5;     { nearly 0.0 - for comparison of reals }
  27.  
  28. Type buffer     = Array [1..maxbuf] Of real;
  29.      trigtable  = Array [0..trigtablesize] Of real;
  30.      contbuffer = Array [1..maxbuf] Of boolean;
  31.      ptypes = (none, mercator, ortho, lambert, azinorth, azisouth);
  32.  
  33. Var projtype : ptypes;              { projection type }
  34.     fasttrig : boolean;             { do we have trig tables? }
  35.     midlon, midlat : real;          { median point for ortho/lambert }
  36.     lonmin, lonmax, latmin, latmax : real; { window borders }
  37.  
  38.  
  39. Function  rmin(x, y : real) : real;          { minimum of two reals }
  40. Function  rmax(x, y : real) : real;          { maximum of two relas }
  41. Function  isignum(i : integer) : integer;    { sign of input (integer) }
  42. Function  rsignum(x : real) : integer;       { sign of input (real) }
  43. Function  scalex(x : real) : integer;        { scale in x direction }
  44. Function  scaley(y : real) : integer;        { scale in y direction }
  45. Function  mercproj(lat:real) : real;         { Mercator projection }
  46. Function  invmercproj(y:real) : real;        { inverse of same }
  47. Procedure orthoproj(lon,lat:real; Var xg,yg,zg:real);{ orthogr. projection }
  48. Procedure orthorot(xg,yg,zg:real; Var x,y,z:real);   { globe rotation }
  49. Procedure lambproj(lon,merclat:real; Var x,y:real);  { Lambert projection }
  50. Procedure aziproj(lon,lat:real; Var x,y:real);       { azimuthal projection }
  51. Procedure project(lon,lat:real; Var x,y:real);       { general projection }
  52. Procedure invproject(x,y:real; Var lon,lat:real);    { inverse of them all }
  53. Procedure checkwindow;                       { correct window, if necessary }
  54. Procedure adaptscale;                        { expand limits to fill screen }
  55. Procedure getpoint(Var x,y:real;             { get coord. via crosshairs }
  56.                    Var finish,firstin:boolean; endgetpoint:boolean);
  57. Procedure drawgrid(gridlon,gridlat:real; grcolour:integer;
  58.                    smile:boolean);           { draw grid on globe }
  59. Procedure mapborder(brcolour:integer; smile:boolean);{ determine&draw borders }
  60. Procedure setprojtype(pt:ptypes);            { set projection parameters }
  61. Procedure readtrigs(trigname : string);      { try to read trig tables }
  62. Procedure dispotrigs;                        { dispose of trig tables }
  63.  
  64.  
  65. Implementation
  66.  
  67. Var
  68.  xlo, xhi, xppu, ylo, yhi, yppu : real;     { scaling parameters }
  69.  zproj : real;                              { z-coordinate of transformation }
  70.  raddeg, pihalf, pifourth, pidegree : real; { Pi-related constants }
  71.  p1,p2,p3,p4,p5,p6,p7,p8,p9 : real;         { projection matrix }
  72.  q1,q2,q3,q4,q5,q6,q7,q8,q9 : real;         { inverse proj. matrix }
  73.  mercptr, sinptr : ^trigtable;            { precomputed Mercator / sine values }
  74.  
  75. Function rmin(x, y : real) : real;
  76. { minimum of arguments                                                       }
  77. Begin                                                                 { rmin }
  78.   If x < y Then rmin := x Else rmin := y;
  79. End;                                                                  { rmin }
  80.  
  81. Function rmax(x, y : real) : real;
  82. { maximum of arguments                                                       }
  83. Begin                                                                 { rmax }
  84.   If x > y Then rmax := x Else rmax := y;
  85. End;                                                                  { rmax }
  86.  
  87. Procedure hilo(Var x, y : buffer; Var c : contbuffer; n : integer);
  88. { update x and y limits of curve buffer                                      }
  89.   Var i : integer;
  90. Begin                                                                 { hilo }
  91.   For i := 1 to n Do
  92.   Begin
  93.     If c[i] Then
  94.     Begin
  95.       If x[i] < xlo Then xlo := x[i];
  96.       If x[i] > xhi Then xhi := x[i];
  97.       If y[i] < ylo Then ylo := y[i];
  98.       If y[i] > yhi Then yhi := y[i];
  99.     End;
  100.   End;
  101. End;                                                                  { hilo }
  102.  
  103. Function isignum(i : integer) : integer;
  104. { -1, if negative; +1, if positive; 0, if 0                                  }
  105. Begin                                                              { isignum }
  106.   If i = 0 Then isignum := 0
  107.            Else If i > 0 Then isignum := 1 Else isignum := -1;
  108. End;                                                               { isignum }
  109.  
  110. Function rsignum(x : real) : integer;
  111. { -1, if negative; +1, if positive; 0, if 0                                  }
  112. Begin                                                              { rsignum }
  113.   If Abs(x) < epsilon Then rsignum := 0
  114.                       Else If x > 0.0 Then rsignum := 1 Else rsignum := -1;
  115. End;                                                               { rsignum }
  116.  
  117. Function scalex(x : real) : integer;
  118. { scale x value to screen dimension                                          }
  119. Begin                                                               { scalex }
  120.   scalex := Trunc(xppu * (x - xlo));
  121. End;                                                                { scalex }
  122.  
  123. Function scaley(y : real) : integer;
  124. { scale y value to screen dimension                                          }
  125. Begin                                                               { scaley }
  126.   scaley := Trunc(yppu * (yhi - y));
  127. End;                                                                { scaley }
  128.  
  129. Function tan(x : real) : real;
  130. { Niklaus forgot it, tsk, tsk                                                }
  131. Begin                                                                  { tan }
  132.   tan := Sin(x) / Cos(x);
  133. End;                                                                   { tan }
  134.  
  135. Function  arcsin(x : real) : real;
  136. Begin                                                               { arcsin }
  137.   arcsin := ArcTan(x / Sqrt(1.0 - Sqr(x)));
  138. End;                                                                { arcsin }
  139.  
  140. Function fastlntan(x : real) : real;
  141. { get Ln(Tan(..)) out of precomputed table                                   }
  142. Begin                                                            { fastlntan }
  143.   If x >= 0.0 Then fastlntan :=  mercptr^[Round(x*100.0)]
  144.               Else fastlntan := -mercptr^[Round(-x*100.0)];
  145. End;                                                             { fastlntan }
  146.  
  147. Function fastsin(x : real) : real;
  148. { get Sin(x) out of precomputed table                                        }
  149.   Var ix : integer;
  150. Begin                                                              { fastsin }
  151.   While x >=  180.0 Do x := x - 360.0;
  152.   While x <= -180.0 Do x := x + 360.0;
  153.   ix := Round(x*100.0);
  154.   If ix >= 0 Then If ix <=  9000 Then fastsin :=  sinptr^[      ix]
  155.                                  Else fastsin :=  sinptr^[18000-ix]
  156.              Else If ix >= -9000 Then fastsin := -sinptr^[     -ix]
  157.                                  Else fastsin := -sinptr^[18000+ix];
  158. End;                                                               { fastsin }
  159.  
  160. Function fastcos(x : real) : real;
  161. { get Cos(x) out of precomputed table                                        }
  162.   Var ix : integer;
  163. Begin                                                              { fastcos }
  164.   While x >=  180.0 Do x := x - 360.0;
  165.   While x <= -180.0 Do x := x + 360.0;
  166.   ix := Round(x*100.0);
  167.   If ix >= 0 Then If ix <=  9000 Then fastcos :=  sinptr^[ 9000-ix]
  168.                                  Else fastcos := -sinptr^[-9000+ix]
  169.              Else If ix >= -9000 Then fastcos :=  sinptr^[ 9000+ix]
  170.                                  Else fastcos := -sinptr^[-9000-ix];
  171. End;                                                               { fastcos }
  172.  
  173. Function mercproj(lat : real) : real;
  174. { compute Mercator projection of latitude                                    }
  175. Begin                                                             { mercproj }
  176.   If fasttrig Then mercproj := fastlntan(lat)
  177.   Else
  178.   Begin
  179.     { leave it alone if too close to pole }
  180.     If (Abs(lat) > 85.0) Then mercproj := lat
  181.                          Else mercproj := Ln(Tan(pidegree*lat + pifourth));
  182.   End;
  183. End;                                                              { mercproj }
  184.  
  185. Function invmercproj(y : real) : real;
  186. { compute inverse of Mercator projection of latitude                         }
  187. Begin                                                          { invmercproj }
  188.   If (Abs(y) > 85.0) Then invmercproj := y
  189.                      Else invmercproj := (ArcTan(Exp(y)) - pifourth)/pidegree;
  190. End;                                                           { invmercproj }
  191.  
  192. Procedure orthoproj(lon, lat : real; Var xg, yg, zg : real);
  193. { calculate R3 coordinates of given point                                    }
  194.   Var coslat : real;
  195. Begin                                                            { orthoproj }
  196.   If fasttrig Then
  197.   Begin
  198.     coslat := fastcos(lat);
  199.     xg := fastcos(lon)*coslat;
  200.     yg := fastsin(lon)*coslat;
  201.     zg := fastsin(lat);
  202.   End
  203.   Else
  204.   Begin
  205.     lon := lon*raddeg;
  206.     lat := lat*raddeg;
  207.     coslat := Cos(lat);
  208.     xg := Cos(lon)*coslat;
  209.     yg := Sin(lon)*coslat;
  210.     zg := Sin(lat);
  211.   End;
  212. End;                                                             { orthoproj }
  213.  
  214. Procedure invorthoproj(xg, yg, zg : real; Var lon, lat : real);
  215. { calculate longitude and latitude                                           }
  216.   Var r : real;
  217. Begin                                                         { invorthoproj }
  218.   r := Sqr(xg) + Sqr(yg);
  219.   r := SqRt(rmax(Sqr(xg)+Sqr(yg),0.0));
  220.   If r > epsilon Then lat := ArcTan(zg/r)/raddeg
  221.                  Else lat := rsignum(zg) * 90.0;
  222.   If Abs(xg) > epsilon Then
  223.   Begin
  224.     lon :=  ArcTan(yg/xg)/raddeg;
  225.     If xg < 0.0 Then If yg >= 0.0 Then lon := lon + 180.0
  226.                                   Else lon := lon - 180.0;
  227.   End Else lon := rsignum(yg) * 90.0;
  228. End;                                                          { invorthoproj }
  229.  
  230. Procedure orthorot(xg, yg, zg : real; Var x, y, z : real);
  231. { rotate globe                                                               }
  232. Begin                                                             { orthorot }
  233.   x := p4*xg + p5*yg;
  234.   y := p7*xg + p8*yg + p9*zg;
  235.   z := p1*xg + p2*yg + p3*zg;
  236. End;                                                              { orthorot }
  237.  
  238. Procedure lambproj(lon, merclat : real; Var x, y : real);
  239. { calculate conformal Lambert mapping of point                               }
  240.   Var alpha, dist : real;
  241. Begin                                                             { lambproj }
  242.   dist := p7 * Exp(-p6 * merclat);
  243.   If fasttrig Then
  244.   Begin
  245.     alpha :=    p6 * (lon-midlon);
  246.     x     :=  dist * fastsin(alpha);
  247.     y     := -dist * fastcos(alpha);
  248.   End
  249.   Else
  250.   Begin
  251.     alpha :=    p4 * (lon-midlon);
  252.     x     :=  dist * Sin(alpha);
  253.     y     := -dist * Cos(alpha);
  254.   End;
  255. End;                                                              { lambproj }
  256.  
  257. Procedure invlambproj(x, y : real; Var lon, lat : real);
  258. { calculate inverse of conformal Lambert mapping                             }
  259.   Var alpha, dist : real;
  260. Begin                                                          { invlambproj }
  261.   If (Abs(x) <= epsilon) And (Abs(y) <= epsilon) Then
  262.   Begin
  263.     lat := 90.0;
  264.     lon :=  0.0;
  265.   End Else
  266.   Begin
  267.     If Abs(y) >= epsilon Then
  268.     Begin
  269.       alpha := ArcTan(-x/y);
  270.       dist  := -y / Cos(alpha);
  271.     End Else
  272.     Begin
  273.       dist := Abs(x);
  274.       alpha := rsignum(x) * 90.0;
  275.     End;
  276.     lon := alpha/p4 + midlon;
  277.     lat := (pifourth - ArcTan(Exp(Ln(dist/p7)/p6)))/pidegree;
  278.   End;
  279. End;                                                           { invlambproj }
  280.  
  281. Procedure aziproj(lon, lat : real; Var x, y : real);
  282. { calculate coordinates under area-prserving azimuthal projection            }
  283.   Var dist : real;
  284. Begin                                                              { aziproj }
  285.   If fasttrig Then
  286.   Begin
  287.     dist := 2.0 * fastsin(45.0 - 0.5*lat);
  288.     x    := dist * fastsin(lon);
  289.     y    := -dist * fastcos(lon);
  290.   End Else
  291.   Begin
  292.     dist := 2.0 * Sin(pifourth - lat*pidegree);
  293.     lon  := lon * raddeg;
  294.     x    := dist * Sin(lon);
  295.     y    := -dist * Cos(lon);
  296.   End;
  297. End;                                                               { aziproj }
  298.  
  299. Procedure invaziproj(x, y : real; Var lon, lat : real);
  300. { inverse azimuthal projection                                               }
  301.   Var dist : real;
  302. Begin                                                           { invaziproj }
  303.   If Abs(y) > epsilon Then
  304.   Begin
  305.     lon  := ArcTan(-x/y) / raddeg;
  306.     If y > 0.0 Then If x < 0.0 Then lon := lon - 180.0
  307.                                Else lon := lon + 180.0;
  308.     dist := SqRt(Sqr(x) + Sqr(y));
  309.   End Else
  310.   Begin
  311.     lon  := rsignum(x) * 90.0;
  312.     dist := Abs(x);
  313.   End;
  314.   lat := 90.0 - ArcSin(0.5*dist) / pidegree;
  315. End;                                                            { invaziproj }
  316.  
  317. Procedure project(lon, lat : real; Var x, y : real);
  318. { convert longitude and latitude to rectangular coordinates                  }
  319.   Var xg, yg, zg: real;   { coordinates of point on globe }
  320. Begin                                                              { project }
  321.   Case projtype Of
  322.     none:      Begin
  323.                  x := lon;
  324.                  y := lat;
  325.                End;
  326.     mercator:  Begin
  327.                  x := lon;
  328.                  y := mercproj(lat);
  329.                End;
  330.     ortho:     Begin
  331.                  orthoproj(lon,lat,xg,yg,zg);
  332.                  { rotate the globe }
  333.                  orthorot(xg,yg,zg,x,y,zproj);
  334.                End;
  335.     lambert:   lambproj(lon,mercproj(lat),x,y);
  336.     azinorth:  aziproj(lon,lat,x,y);
  337.     azisouth:  aziproj(180.0-lon,-lat,x,y);
  338.   End; { Case }
  339. End;                                                               { project }
  340.  
  341. Procedure invproject(x, y : real; Var lon, lat : real);
  342. { inverse projection                                                         }
  343.   Var xg, yg, zg : real;            { coordinates on globe }
  344. Begin                                                           { invproject }
  345.   Case projtype Of
  346.     none:      Begin
  347.                  lon := x;
  348.                  lat := y;
  349.                End;
  350.     mercator:  Begin
  351.                  lon := x;
  352.                  lat := invmercproj(y);
  353.                End;
  354.     ortho:     Begin
  355.                  zproj := SqRt(rmax(1.0-Sqr(x)-Sqr(y),0.0));
  356.                  { rotate the globe }
  357.                  xg := q1*zproj + q2*x + q3*y;
  358.                  yg := q4*zproj + q5*x + q6*y;
  359.                  zg := q7*zproj +        q9*y;
  360.                  invorthoproj(xg,yg,zg,lon,lat);
  361.                End;
  362.     lambert:   invlambproj(x,y,lon,lat);
  363.     azinorth:  invaziproj(x,y,lon,lat);
  364.     azisouth:  Begin
  365.                  invaziproj(x,y,lon,lat);
  366.                  If lon >= 0.0 Then lon :=  180.0 - lon
  367.                                Else lon := -180.0 - lon;
  368.                  lat := -lat;
  369.                End;
  370.   End; {case}
  371. End;                                                            { invproject }
  372.  
  373. Procedure checkwindow;
  374. { adjust window to avoid log(0) and other nasty things                       }
  375.  
  376.   Procedure rswap(Var a, b : real);
  377.   { swap two reals                                                           }
  378.     Var temp : real;
  379.   Begin                                                              { rswap }
  380.     temp := a;
  381.     a := b;
  382.     b := temp;
  383.   End;                                                               { rswap }
  384.  
  385. Begin                                                          { checkwindow }
  386.   If latmin > latmax Then rswap(latmin,latmax);
  387.   If lonmin > lonmax Then rswap(lonmin,lonmax);
  388.   If Abs(latmax-latmin) < epsilon Then
  389.   Begin
  390.     latmax := latmax + 1.0;
  391.     latmin := latmin - 1.0;
  392.   End;
  393.   If Abs(lonmax-lonmin) < epsilon Then
  394.   Begin
  395.     lonmax := lonmax + 1.0;
  396.     lonmin := lonmin - 1.0;
  397.   End;
  398. End;                                                           { checkwindow }
  399.  
  400. Procedure adaptscale;
  401. { expand one of the directions to fill entire screen                         }
  402.   Var mean, scalfac, mlatmin, mlatmax : real;
  403. Begin                                                           { adaptscale }
  404.   Case projtype Of
  405.     none : Begin
  406.              scalfac := xmaxpix / (ymaxpix*aspect);
  407.              If (lonmax-lonmin)/(latmax-latmin) > scalfac Then
  408.              Begin
  409.                { X scale larger than Y: increase Y }
  410.                mean := (latmin + latmax) * 0.5;
  411.                latmin := mean - (lonmax - lonmin) * 0.5 / scalfac;
  412.                latmax := mean + (lonmax - lonmin) * 0.5 / scalfac;
  413.              End Else
  414.              Begin
  415.                { Y scale larger than X: increase X }
  416.                mean := (lonmin + lonmax) * 0.5;
  417.                lonmin := mean - (latmax - latmin) * 0.5 * scalfac;
  418.                lonmax := mean + (latmax - latmin) * 0.5 * scalfac;
  419.              End;
  420.            End;
  421.     mercator : Begin
  422.              checkwindow;
  423.              scalfac := xmaxpix / (ymaxpix*aspect*raddeg);
  424.              mlatmax := mercproj(latmax);
  425.              mlatmin := mercproj(latmin);
  426.              If (lonmax-lonmin)/(mlatmax-mlatmin) > scalfac Then
  427.              Begin
  428.                { X scale larger than Y: increase Y }
  429.                mean := (mlatmin + mlatmax) * 0.5;
  430.                latmin := invmercproj(mean - (lonmax-lonmin)*0.5/scalfac);
  431.                latmax := invmercproj(mean + (lonmax-lonmin)*0.5/scalfac);
  432.                checkwindow;
  433.              End Else
  434.              Begin
  435.                { Y scale larger than X: increase X }
  436.                mean := (lonmin + lonmax) * 0.5;
  437.                lonmin := mean - (mlatmax-mlatmin)*0.5*scalfac;
  438.                lonmax := mean + (mlatmax-mlatmin)*0.5*scalfac;
  439.              End;
  440.            End;
  441.     ortho    : ;  { not applicable }
  442.     lambert  : ;  { not applicable }
  443.     azinorth : ;  { not applicable }
  444.     azisouth : ;  { not applicable }
  445.   End;
  446. End;                                                            { adaptscale }
  447.  
  448. Procedure getlonline(lon : real; Var x, y : buffer;
  449.                                  Var c : contbuffer; Var n : integer);
  450. { create rect. coordinates of line of constant longitude                     }
  451.   Var i : integer;
  452.       lat, step, xx, yy : real;
  453. Begin                                                           { getlonline }
  454.   Case projtype Of
  455.     none, mercator, lambert,
  456.     azinorth, azisouth : Begin
  457.                       project(lon,latmin,x[1],y[1]);
  458.                       project(lon,latmax,x[2],y[2]);
  459.                       c[1] := True; c[2] := True;
  460.                       n := 2;
  461.                     End;
  462.     ortho         : Begin
  463.                       step := (latmax-latmin) / Pred(maxbuf);
  464.                       lat := latmin;
  465.                       For i := 1 To maxbuf Do
  466.                       Begin
  467.                         project(lon,lat,xx,yy);
  468.                         x[i] := xx;
  469.                         y[i] := yy;
  470.                         c[i] := zproj > -epsilon;
  471.                         lat := lat + step;
  472.                       End;
  473.                       n := maxbuf;
  474.                     End;
  475.   End; {case}
  476. End;                                                            { getlonline }
  477.  
  478. Procedure getlatline(lat : real; Var x, y : buffer;
  479.                                  Var c : contbuffer; Var n : integer);
  480. { create rect. coordinates of line of constant latitude                      }
  481.   Var i: integer;
  482.       lon, step, xx, yy : real;
  483. Begin                                                           { getlatline }
  484.   Case projtype Of
  485.     none, mercator: Begin
  486.                       n := 2;
  487.                       project(lonmin,lat,x[1],y[1]);
  488.                       project(lonmax,lat,x[2],y[2]);
  489.                       c[1] := True; c[2] := True;
  490.                     End;
  491.     ortho, lambert,
  492.     azinorth, azisouth : Begin
  493.                       zproj := 1.0; { dummy for all except ortho }
  494.                       If Abs(lat) >= 89.0 Then
  495.                       Begin
  496.                         step := 360.0  / Pred(maxbuf);
  497.                         lon  := midlon - 180.0;
  498.                       End
  499.                       Else Begin
  500.                         step := (lonmax-lonmin) / Pred(maxbuf);
  501.                         lon  := lonmin;
  502.                       End;
  503.                       For i := 1 to maxbuf Do
  504.                       Begin
  505.                         project(lon,lat,xx,yy);
  506.                         x[i] := xx;
  507.                         y[i] := yy;
  508.                         c[i] := zproj > -epsilon;
  509.                         lon := lon + step;
  510.                       End;
  511.                       n := maxbuf;
  512.                     End;
  513.   End; {case}
  514. End;                                                            { getlatline }
  515.  
  516. Procedure getpoint(Var x, y : real;
  517.                    Var finish, firstin : boolean; endgetpoint : boolean);
  518. { get a coordinate pair given by crosshairs                                  }
  519.   Var i, j : integer;
  520.       c : char;
  521. Begin                                                             { getpoint }
  522.   i := Round(xppu*(x - xlo));
  523.   j := Round(yppu*(yhi - y));
  524.   If firstin Then
  525.   Begin
  526.     vline(i);
  527.     hline(j);
  528.     firstin := False;
  529.   End;
  530.   Repeat
  531.     c := ReadKey;
  532.     If (c = #0) Then c := ReadKey;
  533.     Case c Of
  534.      uparr  : If j > 0 Then Begin hline(j); Dec(j); hline(j); End;
  535.      dnarr  : If j < ymaxpix Then Begin hline(j); Inc(j); hline(j); End;
  536.      lfarr  : If i > 0 Then Begin vline(i); Dec(i); vline(i); End;
  537.      rtarr  : If i < xmaxpix Then Begin vline(i); Inc(i); vline(i); End;
  538.      '8', cuparr : If j > 10 Then
  539.                                Begin hline(j); j := j-10; hline(j); End;
  540.      '2', cdnarr : If j < ymaxpix-10 Then
  541.                                Begin hline(j); j := j+10; hline(j); End;
  542.      '4', clfarr : If i > 10 Then
  543.                                Begin vline(i); i := i-10; vline(i); End;
  544.      '6' ,crtarr : If i < xmaxpix-10 Then
  545.                                Begin vline(i); i := i+10; vline(i); End;
  546.     End;
  547.     finish := c In [' ','q','Q',ctrlc,cr,esc];
  548.   Until finish Or Not Endgetpoint;
  549.   If finish Then Begin
  550.     vline(i); hline(j);
  551.   End;
  552.   x := xlo + i/xppu;
  553.   y := yhi - j/yppu
  554. End;                                                              { getpoint }
  555.  
  556. Procedure curve(Var x, y: buffer; Var c : contbuffer; n: integer);
  557. { draw curve of n points                                                     }
  558.   Var xo, yo, xn, yn, i : integer;
  559. Begin                                                                { curve }
  560.   xo := scalex(x[1]); yo := scaley(y[1]);
  561.   For i := 2 to n Do
  562.   Begin
  563.     xn := scalex(x[i]); yn := scaley(y[i]);
  564.     If c[Pred(i)] And c[i] Then Line(xo,yo,xn,yn);
  565.     xo := xn; yo := yn;
  566.   End;
  567. End;                                                                 { curve }
  568.  
  569. Procedure dotcurve(Var x, y: buffer; Var c : contbuffer; n:integer);
  570. { draw dotted curve of n points                                              }
  571.   Var xo, yo, xn, yn, i : integer;
  572.       dotflag : boolean;
  573. Begin                                                             { dotcurve }
  574.   dotflag := True;
  575.   xo := scalex(x[1]); yo := scaley(y[1]);
  576.   For i := 2 to n Do
  577.   Begin
  578.     xn := scalex(x[i]); yn := scaley(y[i]);
  579.     If c[Pred(i)] And c[i] Then dotline(xo,yo,xn,yn,dotflag);
  580.     xo := xn; yo := yn;
  581.   End;
  582. End;                                                              { dotcurve }
  583.  
  584. Procedure drawgrid(gridlon, gridlat : real; grcolour : integer;
  585.                    smile : boolean);
  586. { display grid according to set intervals                                    }
  587.   Var lon, lat : real;
  588.       x, y : buffer;
  589.       c : contbuffer;
  590.       n : integer;
  591.       dotfl : boolean;
  592. Begin                                                             { drawgrid }
  593.   SetColor(word(Abs(grcolour)));
  594.   dotfl := grcolour < 0;
  595.   lon := gridlon * Trunc(lonmin/gridlon);
  596.   If lon < lonmin Then lon := lon + gridlon;
  597.   While lon < lonmax Do
  598.   Begin
  599.     If smile Then showprogress(1);
  600.     getlonline(lon,x,y,c,n);
  601.     If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  602.     lon := lon + gridlon;
  603.   End;
  604.   lat := gridlat * Trunc(latmin/gridlat);
  605.   If lat < latmin Then lat := lat + gridlat;
  606.   While lat < latmax Do
  607.   Begin
  608.     If smile Then showprogress(1);
  609.     getlatline(lat,x,y,c,n);
  610.     If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  611.     lat := lat + gridlat;
  612.   End;
  613.   SetColor(colourglb);
  614. End;                                                              { drawgrid }
  615.  
  616. Procedure mapborder(brcolour : integer; smile : boolean);
  617. { determine scale and draw border of map area                                }
  618.  
  619.   Var x, y, xa, ya  : buffer;
  620.       c, ca : contbuffer;
  621.       n, na : integer;
  622.  
  623.   Procedure mapextremes;
  624.   { determine x and y extremes                                               }
  625.     Var lon, lat, x1, y1, rfact : real;
  626.         i : integer;
  627.   Begin                                                        { mapextremes }
  628.     xlo := 1e30;
  629.     xhi := -1e30;
  630.     ylo := 1e30;
  631.     yhi := -1e30;
  632.     Case projtype Of
  633.       none, mercator, lambert,
  634.       azinorth, azisouth : Begin
  635.           getlatline(latmin,x,y,c,n);
  636.           getlatline(latmax,xa,ya,ca,na);
  637.         End;
  638.       ortho : Begin
  639.           getlonline(lonmin,x,y,c,n);
  640.           getlonline(lonmax,xa,ya,ca,na);
  641.           hilo(x,y,c,n);
  642.           hilo(xa,ya,ca,na);
  643.           If smile Then showprogress(1);
  644.           getlatline(latmin,x,y,c,n);
  645.           getlatline(latmax,xa,ya,ca,na);
  646.           If smile Then showprogress(1);
  647.           hilo(x,y,c,n);
  648.           hilo(xa,ya,ca,na);
  649.           { calculate horizon }
  650.           rfact := 2.0 / Pred(maxbuf);
  651.           x1 := -1.0;
  652.           For i := 1 To maxbuf Do
  653.           Begin
  654.             y1 := 1.0 - Sqr(x1);
  655.             If y1 > 0.0 Then y1 := SqRt(y1) Else y1 := 0.0;
  656.             invproject(x1,y1,lon,lat);
  657.             x[i] := x1; y[i] := y1;
  658.             c[i] :=  (lon >= lonmin-epsilon) And (lon <= lonmax+epsilon) And
  659.                      (lat >= latmin-epsilon) And (lat <= latmax+epsilon);
  660.             invproject(x1,-y1,lon,lat);
  661.             xa[i] := x1; ya[i] := -y1;
  662.             ca[i] := (lon >= lonmin-epsilon) And (lon <= lonmax+epsilon) And
  663.                      (lat >= latmin-epsilon) And (lat <= latmax+epsilon);
  664.             x1 := x1 + rfact;
  665.           End;
  666.           n := maxbuf; na := maxbuf;
  667.           If smile Then showprogress(1);
  668.         End;
  669.     End;
  670.     hilo(x,y,c,n);
  671.     hilo(xa,ya,ca,na);
  672.     If xhi - xlo < epsilon Then xhi := xhi + 1.0;
  673.     If yhi - ylo < epsilon Then yhi := yhi + 1.0;
  674.     If smile Then showprogress(1);
  675.   End;                                                         { mapextremes }
  676.  
  677.   Procedure drawborder;
  678.   { draw border                                                              }
  679.     Var dotfl : boolean;
  680.   Begin                                                         { drawborder }
  681.     SetColor(word(Abs(brcolour)));
  682.     dotfl := brcolour < 0;
  683.     If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  684.     If dotfl Then dotcurve(xa,ya,ca,na) Else curve(xa,ya,ca,na);
  685.     If smile Then showprogress(1);
  686.     If projtype In [ortho,azinorth,azisouth] Then
  687.     Begin
  688.       If lonmax - lonmin < 360.0 - epsilon Then
  689.       Begin
  690.         getlonline(lonmin,x,y,c,n);
  691.         If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  692.         getlonline(lonmax,x,y,c,n);
  693.         If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  694.       End;
  695.       If smile Then showprogress(1);
  696.       getlatline(latmin,x,y,c,n);
  697.       getlatline(latmax,xa,ya,ca,na);
  698.       If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  699.       If dotfl Then dotcurve(xa,ya,ca,na) Else curve(xa,ya,ca,na);
  700.     End
  701.     Else Begin
  702.       getlonline(lonmin,x,y,c,n);
  703.       If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  704.       getlonline(lonmax,x,y,c,n);
  705.       If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
  706.     End;
  707.     SetColor(colourglb);
  708.   End;                                                          { drawborder }
  709.  
  710. Begin                                                            { mapborder }
  711.   mapextremes;
  712.   { set scale parameters (pixels per unit), ensure reasonable aspect ratio }
  713.   xppu := xmaxpix / (xhi-xlo);
  714.   yppu := ymaxpix / (yhi-ylo);
  715.   Case projtype Of
  716.     none, ortho, lambert,
  717.     azinorth, azisouth :  { equal scale in x and y }
  718.              Begin
  719.                xppu := rmin(xppu,yppu*aspect);
  720.                yppu := rmin(yppu,xppu/aspect);
  721.              End;
  722.     mercator : { ppu ratio y/x must be pi/180 }
  723.              Begin
  724.                xppu := rmin(xppu,yppu*aspect*raddeg);
  725.                yppu := rmin(yppu,xppu/(aspect*raddeg));
  726.              End;
  727.   End; {case}
  728.   drawborder;
  729. End;                                                             { mapborder }
  730.  
  731. Procedure setprojtype(pt : ptypes);
  732. { set projection type; calculate parameters                                  }
  733.  
  734.   Procedure lambmatrix;
  735.   { calculate matrix elements for Lambert projection and its inverse         }
  736.   Begin                                                         { lambmatrix }
  737.     p7 := pihalf - midlat*raddeg;
  738.     p6 := Cos(p7);
  739.     p4 := p6 * raddeg;
  740.     p7 := Tan(p7) / Exp(p6 * Ln(Tan(p7*0.5)));
  741.   End;                                                          { lambmatrix }
  742.  
  743.   Procedure orthomatrix;
  744.   { calculate matrix elements for orthogonal projection and its inverse      }
  745.   Begin                                                        { orthomatrix }
  746.     p1 := Cos(midlon*raddeg)*Cos(midlat*raddeg);
  747.     p2 := Sin(midlon*raddeg)*Cos(midlat*raddeg);
  748.     p3 := Sin(midlat*raddeg);
  749.     p4 := -Sin(midlon*raddeg);
  750.     p5 := Cos(midlon*raddeg);
  751.     p6 := 0.0;
  752.     p7 := -Cos(midlon*raddeg)*Sin(midlat*raddeg);
  753.     p8 := -Sin(midlon*raddeg)*Sin(midlat*raddeg);
  754.     p9 := Cos(midlat*raddeg);
  755.     { inverse matrix is transpose}
  756.     q1 := p1;
  757.     q2 := p4;
  758.     q3 := p7;
  759.     q4 := p2;
  760.     q5 := p5;
  761.     q6 := p8;
  762.     q7 := p3;
  763.     q8 := 0.0;
  764.     q9 := p9;
  765.   End;                                                         { orthomatrix }
  766.  
  767. Begin                                                             { projtype }
  768.   projtype := pt;
  769.   Case projtype Of
  770.     none     : ;
  771.     mercator : ;
  772.     ortho    : orthomatrix;
  773.     lambert  : lambmatrix;
  774.     azinorth, azisouth : Begin
  775.                  If latmax >= -latmin Then projtype := azinorth
  776.                                       Else projtype := azisouth;
  777.                  midlon := 0.0; midlat := 0.0;
  778.                End;
  779.   End;
  780. End;                                                              { projtype }
  781.  
  782. Procedure readtrigs(trigname : string);
  783. { read trig tables into memory, if possible                                  }
  784.   Var trigf : File Of trigtable;
  785. Begin                                                            { readtrigs }
  786.   New(mercptr);
  787.   New(sinptr);
  788.   fasttrig := (mercptr <> Nil) And (sinptr <> Nil);
  789.   If (mercptr <> Nil) And (sinptr = Nil) Then Dispose(mercptr);
  790.   If fasttrig Then
  791.   Begin
  792.     Assign(trigf,trigname);
  793.     {$I- } Reset(trigf);
  794.     read(trigf,mercptr^);
  795.     read(trigf,sinptr^);
  796.     Close(trigf);
  797.     {$I+ }
  798.     fasttrig := IOResult = 0;
  799.   End;
  800. End;                                                             { readtrigs }
  801.  
  802. Procedure dispotrigs;
  803. { dispose of trig tables in memory                                           }
  804. Begin                                                           { dispotrigs }
  805.   If fasttrig Then
  806.   Begin
  807.     Dispose(mercptr);
  808.     Dispose(sinptr);
  809.     fasttrig := False;
  810.   End;
  811. End;                                                            { dispotrigs }
  812.  
  813. Begin
  814.   fasttrig := False;
  815.   raddeg   := Pi / 180.0;
  816.   pihalf   := Pi / 2.0;
  817.   pifourth := Pi / 4.0;
  818.   pidegree := Pi / 360.0;
  819.   projtype := mercator;
  820.   midlon   := 0.0;
  821.   midlat   := 0.0;
  822.   lonmin   := -200.0;
  823.   lonmax   :=  200.0;
  824.   latmin   :=  -80.0;
  825.   latmax   :=   80.0;
  826.   xppu     := 1.0;
  827.   yppu     := 1.0;
  828. End.
  829.