home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d3xx / d322 / gwin.lha / Gwin / examples / geomap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-02-27  |  12.3 KB  |  594 lines
  1. #include "gwin.user.h"
  2. #include <math.h>
  3. double myatan2();
  4. float t0[3],t1[3],t2[3],ti0[3],ti1[3],ti2[3];
  5. double pi,pi1, pi2, pi3;
  6. float xmerc();
  7. float lats[500],longs[500];
  8. int ortho;
  9. int merc;
  10. int setmercator(),setortho(),plot_grid(),interrupt_display();
  11. int new_center(),redisplay(),setortho2();
  12. float x,y;
  13. int event;
  14. char key;
  15. int interrupt = 0;
  16. int display_underway = 0;
  17. float lat_cen,long_cen;
  18. int ortho2;
  19.  
  20. main()
  21. {
  22. float xlat,xlong;
  23. char line[30];
  24.  
  25.    pi1 = 3.1415926535897943/180.;
  26.    pi2 = 180./3.1415926535897943;
  27.    pi3 = 3.1415926535897943/2.0;
  28.  
  29.    ortho = 1;
  30.    ortho2 = 0;
  31.    merc = 0;
  32.  
  33.    printf("\n\nEnter lat_cen, long_cen:  ");
  34.    scanf("%f %f",&lat_cen,&long_cen);
  35.  
  36.    ustart("high2",0.,640.,0.,400.);
  37.    if(merc ) uwindo(-180.,180.,-180.,180.);
  38.    if(ortho || ortho2) uwindo(-250.,250.,-160.,160.);
  39.  
  40.    setup_menus();
  41.  
  42.    make_rotation_matrices(lat_cen,long_cen);
  43.  
  44.    make_geo_grid();
  45.  
  46.    plot_grid();
  47.  
  48.    while(1 == 1) {
  49.       ugrina(&x,&y,&event,&key);
  50.       latlong(x,y,&xlat,&xlong);
  51.  
  52.       sprintf(line,"%6.2f  %6.2f",xlat,xlong);
  53.       if (ortho || ortho2){
  54.          uprint(100.,150.,line);
  55.       }else{
  56.          uprint(100.,170.,line);
  57.       }
  58.  
  59.    }
  60. }
  61.  
  62. latlong(x1,y1,xlat,xlong)
  63. float x1,y1,*xlat,*xlong;
  64. {
  65. float xxx,yyy,zzz;
  66. double xtemp,ytemp,ztemp;
  67. float temp1,temp2;
  68. float ph1,th1;
  69. char line[30];
  70.  
  71.    if(ortho || ortho2){
  72.       yyy = x1/150.;
  73.       zzz = y1/150.;
  74.  
  75.       temp1 = yyy*yyy;
  76.       temp2 = zzz*zzz;
  77.       if(temp1 + temp2 > 1.0)return(1);
  78.  
  79.       xxx = sqrt((double)(1 - temp1 - temp2));
  80.  
  81.       xtemp = (double)(ti0[0] * xxx + ti0[1] * yyy + ti0[2] * zzz);
  82.       ytemp = (double)(ti1[0] * xxx + ti1[1] * yyy + ti1[2] * zzz);
  83.       ztemp = (double)(ti2[0] * xxx + ti2[1] * yyy + ti2[2] * zzz);
  84.  
  85.       if(fabs(xtemp) < 1e-15) xtemp = 1e-15;
  86.       *xlong = pi2*myatan2(ytemp,xtemp);
  87.  
  88.       if(ztemp >  1.0) ztemp =  1.0;
  89.       if(ztemp < -1.0) ztemp = -1.0;
  90.       *xlat = pi2*(pi3 - acos(ztemp));
  91.  
  92.    }
  93.  
  94.    if(merc){
  95.  
  96.       ph1 = pi1 * y1;
  97.       ph1 = 2.0 * atan(exp(ph1)) - pi3;
  98.       ph1 = pi3 - ph1;
  99.  
  100.       th1 = pi1 * x1;
  101.  
  102.       xxx = sin((double)ph1) * cos((double)th1);
  103.       yyy = sin((double)ph1) * sin((double)th1);
  104.       zzz = cos((double)ph1);
  105.  
  106.       xtemp = (double)(ti0[0] * xxx + ti0[1] * yyy + ti0[2] * zzz);
  107.       ytemp = (double)(ti1[0] * xxx + ti1[1] * yyy + ti1[2] * zzz);
  108.       ztemp = (double)(ti2[0] * xxx + ti2[1] * yyy + ti2[2] * zzz);
  109.  
  110.       if(fabs(xtemp) < 1e-15) xtemp = 1e-15;
  111.       *xlong = pi2*myatan2(ytemp,xtemp);
  112.  
  113.       if(ztemp >  1.0) ztemp =  1.0;
  114.       if(ztemp < -1.0) ztemp = -1.0;
  115.       *xlat = pi2*(pi3 - acos(ztemp));
  116.  
  117.  
  118.    }
  119. }
  120.  
  121. make_rotation_matrices(lat_cen,long_cen)
  122. float lat_cen,long_cen;
  123. {
  124.  
  125. /* Rotation matrices to rotate lat/long to new lat/long */
  126. /* so that (lat_cen,long_cen) becomes (0,0)             */
  127.  
  128. float sin_theta, cos_theta, theta_rot, sin_phi, cos_phi, phi_rot;
  129.  
  130. theta_rot = long_cen * pi1;
  131. phi_rot   = lat_cen  * pi1;
  132.  
  133. sin_theta = sin(theta_rot);
  134. cos_theta = cos(theta_rot);
  135. sin_phi   = sin(phi_rot);
  136. cos_phi   = cos(phi_rot);
  137.  
  138. t0[0] = cos_theta * cos_phi;
  139. t0[1] = sin_theta * cos_phi;
  140. t0[2] = sin_phi;
  141. t1[0] = -sin_theta;
  142. t1[1] = cos_theta;
  143. t1[2] = 0.0;
  144. t2[0] = -cos_theta * sin_phi;
  145. t2[1] = -sin_theta * sin_phi;
  146. t2[2] = cos_phi;
  147.  
  148. ti0[0] = cos_theta * cos_phi;
  149. ti0[1] = -sin_theta;
  150. ti0[2] = -cos_theta * sin_phi;
  151. ti1[0] = sin_theta * cos_phi;
  152. ti1[1] = cos_theta;
  153. ti1[2] = -sin_theta * sin_phi;
  154. ti2[0] = sin_phi;
  155. ti2[1] = 0.0;
  156. ti2[2] = cos_phi;
  157.  
  158. }
  159.  
  160. mercator(long_in,lat_in,x_out,y_out)
  161. float lat_in,long_in,*y_out,*x_out;
  162. {
  163. float xla,xlo,x[3],xr[3],sth,cth,sph,cph;
  164. double th1,ph1;
  165.  
  166.    xla = lat_in;
  167.    xlo = long_in;
  168.    th1 = pi1 * xlo;
  169.    ph1 = pi3 - xla * pi1;
  170.  
  171.    sth = sin(th1);
  172.    cth = cos(th1);
  173.    sph = sin(ph1);
  174.    cph = cos(ph1);
  175.  
  176.    x[0] = sph * cth;
  177.    x[1] = sph * sth;
  178.    x[2] = cph;
  179.  
  180.    xr[0] = t0[0] * x[0] + t0[1] * x[1] + t0[2] * x[2];
  181.    xr[1] = t1[0] * x[0] + t1[1] * x[1] + t1[2] * x[2];
  182.    xr[2] = t2[0] * x[0] + t2[1] * x[1] + t2[2] * x[2];
  183.  
  184.    *y_out = xmerc( pi2 * asin(xr[2]));
  185.    *x_out = pi2 * myatan2(xr[1],xr[0]);
  186. }
  187.  
  188. float xmerc(xlat)
  189. float(xlat);
  190. {
  191. double a1;
  192.  
  193.    if (fabs(xlat) < 89.9999){
  194.       if(xlat >= 0.0){
  195.          a1 = pi1 * (45.0 + .5 * xlat);
  196.          return((float) pi2 * log(tan(a1)));
  197.       } else {
  198.          a1 = -pi1 * (-45.0 + .5 * xlat);
  199.          return((float) -pi2 * log(tan(a1)));
  200.       }
  201.    } else {
  202.  
  203.       if(xlat < 0) return(-750.);
  204.       return(750.);
  205.    }
  206. }
  207.  
  208. orthographic(long_in,lat_in,x_out,y_out,frontside)
  209. float lat_in,long_in,*y_out,*x_out;
  210. int *frontside;
  211. {
  212. float xla,xlo,x[3],xr[3],sth,cth,sph,cph;
  213. double th1,ph1;
  214.  
  215.    xla = lat_in;
  216.    xlo = long_in;
  217.    th1 = pi1 * xlo;
  218.    ph1 = pi3 - xla * pi1;
  219.  
  220.    sth = sin(th1);
  221.    cth = cos(th1);
  222.    sph = sin(ph1);
  223.    cph = cos(ph1);
  224.  
  225.    x[0] = sph * cth;
  226.    x[1] = sph * sth;
  227.    x[2] = cph;
  228.  
  229.    xr[0] = t0[0] * x[0] + t0[1] * x[1] + t0[2] * x[2];
  230.    xr[1] = t1[0] * x[0] + t1[1] * x[1] + t1[2] * x[2];
  231.    xr[2] = t2[0] * x[0] + t2[1] * x[1] + t2[2] * x[2];
  232.  
  233.    *y_out = 150.0 * xr[2];
  234.    *x_out = 150.0 * xr[1];
  235.    *frontside = 0;
  236.    if(xr[0]>0.0)*frontside = 1;
  237.  
  238. }
  239.  
  240. make_geo_grid()
  241. {
  242. int i;
  243. float x;
  244.  
  245.    i = 0;
  246.  
  247.    for(x = -90.0;x < 91.0 ;x += 10.0){
  248.       lats[i++] = x;
  249.    }
  250.    lats[0] = -89.999;
  251.    lats[i-1] = 89.999;
  252.    i = 0;
  253.  
  254.    for(x = -180.0; x < 181.0; x += 10.0){
  255.       longs[i++] = x;
  256.    }
  257.    longs[0] = -179.999;
  258.    longs[i-1] = 179.999;
  259. }
  260.  
  261. plot_grid()
  262. {
  263.    if(merc){
  264.       mercplot();
  265.       return(0);
  266.    }
  267.    if(ortho){
  268.       orthoplot();
  269.       return(0);
  270.    }
  271.    if(ortho2){
  272.       ortho2plot();
  273.       return(0);
  274.    }
  275. }
  276.  
  277. mercplot()
  278. {
  279. int i,j;
  280. float x,y,xold;
  281. int event,key;
  282.  
  283.    display_underway = 1;
  284.    interrupt = 0;
  285.  
  286.    uerase();
  287.  
  288.    upset("colo",1.0);
  289.  
  290.    for(i=0;i<19;i++){
  291.       mercator(longs[0],lats[i],&x,&y);
  292.       umove(x,y);
  293.       xold = x;
  294.       for(j=1;j<37;j++){
  295.          ugrina(&x,&y,&event,&key);
  296.          if(interrupt) goto interrupted;
  297.          mercator(longs[j],lats[i],&x,&y);
  298.          if(fabs(xold-x) < 100.0)udraw(x,y);
  299.          umove(x,y);
  300.          xold = x;
  301.       }
  302.    }
  303.    for(j=0;j<37;j++){
  304.       mercator(longs[j],lats[0],&x,&y);
  305.       umove(x,y);
  306.       xold = x;
  307.       for(i=1;i<19;i++){
  308.          ugrina(&x,&y,&event,&key);
  309.          if(interrupt) goto interrupted;
  310.          mercator(longs[j],lats[i],&x,&y);
  311.          if(fabs(xold-x) < 100.0)udraw(x,y);
  312.          umove(x,y);
  313.          xold = x;
  314.       }
  315.    }
  316. interrupted:
  317.    interrupt = 0;
  318.    display_underway = 0;
  319. }
  320.  
  321. orthoplot()
  322. {
  323. int i,j;
  324. float x,y,xold,yold;
  325. int event,key;
  326. int frontside;
  327.  
  328.    display_underway = 1;
  329.    interrupt = 0;
  330.  
  331.    uerase();
  332.  
  333.    upset("colo",1.0);
  334.  
  335.    for(i=0;i<19;i++){
  336.       orthographic(longs[0],lats[i],&x,&y,&frontside);
  337.       umove(x,y);
  338.       xold = x;
  339.       yold = y;
  340.       for(j=1;j<37;j++){
  341.          ugrina(&x,&y,&event,&key);
  342.          if(interrupt) goto interrupted;
  343.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  344.          if(fabs(xold-x) < 100.0 && frontside){
  345.             umove(xold,yold);
  346.             udraw(x,y);
  347.          }else{
  348.             umove(x,y);
  349.          }
  350.          xold = x;
  351.          yold = y;
  352.       }
  353.    }
  354.    for(j=0;j<37;j++){
  355.       orthographic(longs[j],lats[0],&x,&y,&frontside);
  356.       umove(x,y);
  357.       xold = x;
  358.       yold = y;
  359.       for(i=1;i<19;i++){
  360.          ugrina(&x,&y,&event,&key);
  361.          if(interrupt) goto interrupted;
  362.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  363.          if(fabs(xold-x) < 100.0 && frontside){
  364.             umove(xold,yold);
  365.             udraw(x,y);
  366.          }else{
  367.             umove(x,y);
  368.          }
  369.          xold = x;
  370.          yold = y;
  371.       }
  372.    }
  373. interrupted:
  374.    interrupt = 0;
  375.    display_underway = 0;
  376. }
  377.  
  378. ortho2plot()
  379. {
  380. int i,j;
  381. float x,y,xold,yold;
  382. int event,key;
  383. int frontside;
  384.  
  385.    frontside = 1;
  386.    display_underway = 1;
  387.    interrupt = 0;
  388.  
  389.    uerase();
  390.  
  391. /* plot backsides first */
  392.  
  393.    upset("colo",10.0);
  394.  
  395.    for(i=0;i<19;i++){
  396.       orthographic(longs[0],lats[i],&x,&y,&frontside);
  397.       umove(x,y);
  398.       xold = x;
  399.       yold = y;
  400.       for(j=1;j<37;j++){
  401.          ugrina(&x,&y,&event,&key);
  402.          if(interrupt) goto interrupted;
  403.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  404.          if(fabs(xold-x) < 100.0 && !frontside){
  405.             umove(xold,yold);
  406.             udraw(x,y);
  407.          }else{
  408.             umove(x,y);
  409.          }
  410.          xold = x;
  411.          yold = y;
  412.       }
  413.    }
  414.    for(j=0;j<37;j++){
  415.       orthographic(longs[j],lats[0],&x,&y,&frontside);
  416.       umove(x,y);
  417.       xold = x;
  418.       yold = y;
  419.       for(i=1;i<19;i++){
  420.          ugrina(&x,&y,&event,&key);
  421.          if(interrupt) goto interrupted;
  422.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  423.          if(fabs(xold-x) < 100.0 && !frontside){
  424.             umove(xold,yold);
  425.             udraw(x,y);
  426.          }else{
  427.             umove(x,y);
  428.          }
  429.          xold = x;
  430.          yold = y;
  431.       }
  432.    }
  433.  
  434. /* now plot front sides */
  435.  
  436.    upset("colo",1.0);
  437.  
  438.    for(i=0;i<19;i++){
  439.       orthographic(longs[0],lats[i],&x,&y,&frontside);
  440.       umove(x,y);
  441.       xold = x;
  442.       yold = y;
  443.       for(j=1;j<37;j++){
  444.          ugrina(&x,&y,&event,&key);
  445.          if(interrupt) goto interrupted;
  446.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  447.          if(fabs(xold-x) < 100.0 && frontside){
  448.             umove(xold,yold);
  449.             udraw(x,y);
  450.          }else{
  451.             umove(x,y);
  452.          }
  453.          xold = x;
  454.          yold = y;
  455.       }
  456.    }
  457.    for(j=0;j<37;j++){
  458.       orthographic(longs[j],lats[0],&x,&y,&frontside);
  459.       umove(x,y);
  460.       xold = x;
  461.       yold = y;
  462.       for(i=1;i<19;i++){
  463.          ugrina(&x,&y,&event,&key);
  464.          if(interrupt) goto interrupted;
  465.          orthographic(longs[j],lats[i],&x,&y,&frontside);
  466.          if(fabs(xold-x) < 100.0 && frontside){
  467.             umove(xold,yold);
  468.             udraw(x,y);
  469.          }else{
  470.             umove(x,y);
  471.          }
  472.          xold = x;
  473.          yold = y;
  474.       }
  475.    }
  476. interrupted:
  477.    interrupt = 0;
  478.    display_underway = 0;
  479. }
  480.  
  481. setup_menus()
  482. {
  483.    uamenu(1,0,0,"PROJECTION      ",' ',0,(USHORT)MIDRAWN|MENUENABLED,7);
  484.    uamenu(1,1,0,"   ORTHOGRAPHIC ",'O',6,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  485.           |COMMSEQ|ITEMENABLED|CHECKIT,setortho);
  486.    uamenu(1,2,0,"   MERCATOR     ",'M',5,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  487.           |COMMSEQ|ITEMENABLED,setmercator);
  488.    uamenu(1,3,0,"   ORTHO2       ",'2',3,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  489.           |COMMSEQ|ITEMENABLED,setortho2);
  490.  
  491.    uamenu(2,0,0,"DISPLAY    ",' ',0,(USHORT)MIDRAWN|MENUENABLED,7);
  492.    uamenu(2,1,0,"REDISPLAY  ",'R',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  493.           |COMMSEQ|ITEMENABLED,redisplay);
  494.    uamenu(2,2,0,"INTERRUPT  ",'I',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  495.           |COMMSEQ|ITEMENABLED,interrupt_display);
  496.    uamenu(2,3,0,"NEW CENTER ",'C',0,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  497.           |COMMSEQ|ITEMENABLED,new_center);
  498.  
  499. }
  500.  
  501. redisplay()
  502. {
  503.    plot_grid();
  504.    interrupt=1;  /* cancel previous invocation if any */
  505. }
  506.  
  507. setortho()
  508. {
  509.    uwindo(-250.,250.,-160.,160.);
  510.    ortho = 1;
  511.    merc = 0;
  512.    ortho2 = 0;
  513.    uamenu(1,1,0,"   ORTHOGRAPHIC ",'O',6,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  514.           |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setortho);
  515.    redisplay();
  516. }
  517.  
  518. setmercator()
  519. {
  520.    uwindo(-180.,180.,-180.,180.);
  521.    merc = 1;
  522.    ortho = 0;
  523.    ortho2 = 0;
  524.    uamenu(1,2,0,"   MERCATOR     ",'M',5,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  525.           |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setmercator);
  526.    redisplay();
  527. }
  528.  
  529. setortho2()
  530. {
  531.    uwindo(-250.,250.,-160.,160.);
  532.    ortho2 = 1;
  533.    ortho = 0;
  534.    merc = 0;
  535.    uamenu(1,3,0,"   ORTHO2       ",'2',3,(USHORT)MIDRAWN|ITEMTEXT|HIGHCOMP
  536.           |COMMSEQ|ITEMENABLED|CHECKIT|CHECKED,setortho2);
  537.    redisplay();
  538. }
  539.  
  540. interrupt_display()
  541. {
  542.     interrupt = 1;
  543. }
  544.  
  545. new_center()
  546. {
  547. char coords[255];
  548.  
  549.    uerase();
  550.    uwindo(0.,100.,0.,100.);
  551.    ugetstring(10.,75.,70.,"Enter latitude and longitude:  ",coords);
  552.  
  553.    while (sscanf(coords,"%f %f",&lat_cen,&long_cen) < 2){
  554.       ugetstring(10.,75.,70.,
  555.          "Error, try again.  Enter latitude and longitude:  ",coords);
  556.    }
  557.    make_rotation_matrices(lat_cen,long_cen);
  558.  
  559.    if(merc)setmercator();
  560.    if(ortho)setortho();
  561.    if(ortho2)setortho2();
  562. }
  563.  
  564. double myatan2(y,x)
  565. double x,y;
  566. {
  567. double temp;
  568. double pi = 3.1415926535897943;
  569. double piov2 = 0.5 * 3.1415926535897943;
  570.  
  571.    if(fabs(x) > 0.0){
  572.       temp = atan((double)y/x);
  573.    }else{
  574.       if(y < 0.0){
  575.          temp = -piov2;
  576.       }else{
  577.          temp =  piov2;
  578.          if(fabs(x) < 1e-15 && fabs(y) <1e-15) temp = 0.0;
  579.       }
  580.    }
  581.  
  582.    if(x < 0.0 ){
  583.       if(y >= 0.0) {
  584.          temp =   pi + temp;
  585.       }else{
  586.          temp =  -pi + temp;
  587.       }
  588.    }
  589.    return(temp);
  590. }
  591.  
  592.  
  593.  
  594.