home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / amiglobe099src.lha / Ami-Globe / projection.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-05  |  14.6 KB  |  414 lines
  1. /************************************************************************/
  2. /*                                                                      */
  3. /*      fichier         : projection.c                                  */
  4. /*      projet          : amiglobe                                      */
  5. /*      date création   : 10/09/94                                      */
  6. /*      commentaire     : fonctions de conversion de                    */
  7. /*                        points écran <-> point de la carte            */
  8. /*      révision        : $VER: projection.c 1.008 (30 Sep 1994) */
  9. /*      copyright       : Olivier Collard, Thomas Landspurg             */
  10. /*                                                                      */
  11. /************************************************************************/
  12.  
  13.  
  14. /************************************************************************/
  15. /*      includes                                                        */
  16. /************************************************************************/
  17. #include <math.h>
  18. #include <stdlib.h>
  19. #include <exec/types.h>
  20. #include <stdio.h>
  21. #include <clib/macros.h>
  22. #include <libraries/fastsincos.h>
  23. #include "amiglobe_types.h"
  24. #include "database_types.h"
  25. #include "map_function_protos.h"
  26. #include "3D_types.h"
  27. #include "3D_protos.h"
  28. #include "projection.h"
  29.  
  30. /************************************************************************/
  31. /*      variables externes                                              */
  32. /************************************************************************/
  33. extern PREFERENCE Pref;
  34. extern CLIP clip_proj;
  35. extern CLIP clip_max;
  36. extern struct Library * FastsincosBase;
  37.  
  38. /************************************************************************/
  39. /*      variables privées                                               */
  40. /************************************************************************/
  41. double cos_b2=0;
  42. double sin_b2=0;
  43. int mx=0;
  44. int my=0;
  45. double win_ratio;
  46.  
  47.  
  48. /************************************************************************/
  49. /*      implémentation                                                  */
  50. /************************************************************************/
  51.  
  52. /*************************************************************************/
  53. void conv_clip(void)
  54. {
  55.     switch(Pref.Type_Proj)
  56.     {
  57.         case PROJ_GLOBE:
  58.             clip_proj=Pref.clip_cur;
  59.             /*{
  60.                 int mx=(Pref.clip_cur.minx+Pref.clip_cur.maxx)/2;
  61.                 clip_proj=Pref.clip_cur;
  62.                 if (mx-clip_proj.minx>9000)
  63.                     clip_proj.minx=mx-9000;
  64.                 if (clip_proj.maxx-mx>9000)
  65.                     clip_proj.maxx=mx+9000;
  66.             }*/
  67.             break;
  68.         case PROJ_MERCATOR:
  69.             clip_proj.maxy=
  70.                 sin(Map_Convert_Angle(Pref.clip_cur.maxy))*clip_max.maxy;
  71.             clip_proj.miny=
  72.                 sin(Map_Convert_Angle(Pref.clip_cur.miny))*clip_max.maxy;
  73.             clip_proj.maxx=Pref.clip_cur.maxx;
  74.             clip_proj.minx=Pref.clip_cur.minx;
  75.             break;
  76.         case PROJ_FLAT:
  77.             clip_proj=Pref.clip_cur;
  78.             break;
  79.         default:
  80.             clip_proj=Pref.clip_cur;
  81.     }
  82. }
  83.  
  84. /* ------------------------------ Update_New_Clip ------------------------------
  85.  
  86.  Commentaire:cette fonction recalcule les constantes servant à conv_xy
  87.              et conv_inv_xy à chaque changement du clip_cur
  88.              CETTE FONCTION DOIT ETRE APPELEE A CHAQUE NOUVEAU CLIP_CUR
  89.  
  90. */
  91.  
  92. void
  93. Update_New_Clip(void)
  94. {
  95.     static double b2;
  96.     mx=(Pref.clip_cur.maxx+Pref.clip_cur.minx)/2;
  97.     my=(Pref.clip_cur.maxy+Pref.clip_cur.miny)/2;
  98.     b2=my*PI/18000;
  99.     cos_b2=cos(b2);
  100.     sin_b2=sin(b2);
  101.     win_ratio=Pref.win_width/Pref.win_height;
  102.  
  103. }
  104.  
  105. /****************************************************************************/
  106. /* converti un point écran en un point de la carte*/
  107. void conv_inv_xy(
  108.         int     *px,int *py
  109.         )
  110. {
  111.     double  f_angle;//,f_lat,f_lon;
  112.     int     x,y;
  113.     x=*px;
  114.     y=*py;
  115.  
  116.     if(Pref.Flg_Proj_3D)
  117.     {
  118.         T3D_Convert_Inv(x,y,px,py);
  119.     }
  120.     switch(Pref.Type_Proj)
  121.     {
  122.     case PROJ_GLOBE:
  123.     {
  124.             static float alpha; /* peut être <0 */
  125.  
  126.             static float gamma; /* peut être <0 */
  127.             static float c;     /* [0,PI]      */
  128.             static float a;     /* [0,PI]      */
  129.     x=(x*(Pref.clip_cur.maxx-Pref.clip_cur.minx))/(Pref.win_height*win_ratio)+
  130.             Pref.clip_cur.minx-mx;
  131.     y=(y*(Pref.clip_cur.maxy-Pref.clip_cur.miny))/(Pref.win_height)+
  132.             Pref.clip_cur.miny-my;
  133.             a=acos(-(double)y/(double)clip_max.maxy);
  134.             if (FastsincosBase!=NULL)
  135.             {
  136.                 gamma=asin((double)x/((double)clip_max.maxx*fastsin(a)));
  137.                 if (fabs ((float)y/(float)clip_max.maxy) <=1
  138.                     && fabs( (float)x / ((float)clip_max.maxx*fastcos(a-PI/2)) ) <=1 )
  139.                 {
  140.     /*          if (f_lat>PI/2)
  141.                 {
  142.                     f_lat=PI-f_lat;
  143.                     if (f_lon>0)
  144.                         f_lon-=PI;
  145.                     else
  146.                         f_lon+=PI;
  147.                 }
  148.                 else
  149.                 if (f_lat<PI/2)
  150.                 {
  151.                         f_lat=-PI+f_lat;
  152.                         if (f_lon>0) 
  153.                                 f_lon-=PI;
  154.                         else
  155.                                 f_lon+=PI;
  156.                 }*/     
  157.                     c=acos(fastcos(a)*cos_b2-fastsin(a)*sin_b2*fastcos(gamma));
  158.                     alpha=asin(fastsin(a)*fastsin(gamma)/fastsin(c));
  159.                     *px=alpha*18000/PI+mx;
  160.                     if (*px>18000)
  161.                         *px-=18000;
  162.                     if (*px<-18000)
  163.                         *px+=18000;
  164.                     *py=-9000+c*clip_max.maxy*2/PI; 
  165.                 }
  166.                 else
  167.                 {
  168.                     if (x<0) 
  169.                         *px=-clip_max.maxx;
  170.                     else
  171.                         *px=clip_max.maxx;
  172.                     if (y<0)
  173.                         *py=-clip_max.maxy;
  174.                     else
  175.                         *py=clip_max.maxy;
  176.                 }
  177.             }
  178.             else
  179.             {
  180.                 gamma=asin((double)x/((double)clip_max.maxx*sin(a)));
  181.                 if (fabs ((float)y/(float)clip_max.maxy) <=1
  182.                     && fabs( (float)x / ((float)clip_max.maxx*cos(a-PI/2)) ) <=1 )
  183.                 {
  184.     /*          if (f_lat>PI/2)
  185.                 {
  186.                     f_lat=PI-f_lat;
  187.                     if (f_lon>0)
  188.                         f_lon-=PI;
  189.                     else
  190.                         f_lon+=PI;
  191.                 }
  192.                 else
  193.                 if (f_lat<PI/2)
  194.                 {
  195.                         f_lat=-PI+f_lat;
  196.                         if (f_lon>0) 
  197.                                 f_lon-=PI;
  198.                         else
  199.                                 f_lon+=PI;
  200.                 }*/     
  201.                     c=acos(cos(a)*cos_b2-sin(a)*sin_b2*cos(gamma));
  202.                     alpha=asin(sin(a)*sin(gamma)/sin(c));
  203.                     *px=alpha*18000/PI+mx;
  204.                     *py=-9000+c*clip_max.maxy*2/PI;
  205.                 }
  206.                 else
  207.                 {
  208.                     if (x<0) 
  209.                         *px=-clip_max.maxx;
  210.                     else
  211.                         *px=clip_max.maxx;
  212.                     if (y<0)
  213.                         *py=-clip_max.maxy;
  214.                     else
  215.                         *py=clip_max.maxy;
  216.                 }
  217.             }
  218.     }
  219.     break;
  220.     case PROJ_MERCATOR:
  221.         x=(x*(Pref.clip_cur.maxx-Pref.clip_cur.minx))/(Pref.win_width)+Pref.clip_cur.minx;
  222.         y=(y*(clip_proj.maxy-clip_proj.miny))/(Pref.win_height)+clip_proj.miny;
  223.         if (y>clip_max.maxy)
  224.             y=clip_max.maxy;
  225.         if (y<clip_max.miny)
  226.             y=clip_max.miny;
  227.         f_angle=(double)asin((float)((float)y/(float)clip_max.maxy));
  228.         y=f_angle*(clip_max.maxx-clip_max.minx)/(2*PI);
  229.         *px=x;
  230.         *py=y;          
  231.         break;
  232.  
  233.     case PROJ_FLAT:
  234.     default:
  235.         x=(x*(Pref.clip_cur.maxx-Pref.clip_cur.minx))/(Pref.win_width)+Pref.clip_cur.minx;
  236.         y=(y*(Pref.clip_cur.maxy-Pref.clip_cur.miny))/(Pref.win_height)+Pref.clip_cur.miny;
  237.         *px=x;
  238.         *py=y;
  239.         break;
  240.     }
  241. }
  242.  
  243. /****************************************************************************/
  244. /* converti un point de la carte en point écran*/
  245. void conv_xy(int *px,int *py)
  246. {
  247.     float   f_angle;//,f_lat,f_lon;
  248.     int     x,y;
  249. /*      int clipmaxy,clipminy;*/
  250.  
  251.     x=*px;
  252.     y=*py;
  253.  
  254.     switch (Pref.Type_Proj)
  255.     {
  256.  
  257.     case PROJ_MERCATOR:
  258.             f_angle=Map_Convert_Angle(y);
  259.             y=sin(f_angle)*(clip_max.maxy);
  260. *px=(((x-Pref.clip_cur.minx)*(Pref.win_width))/(Pref.clip_cur.maxx - Pref.clip_cur.minx));
  261. *py=(((y-clip_proj.miny)*(Pref.win_height))/(clip_proj.maxy - clip_proj.miny));
  262.             
  263. /*                      *py=y;
  264.             *px=x;      */
  265.             break;
  266.     case PROJ_GLOBE:
  267.         {
  268.             static float alpha; /* peut être <0 */
  269.  
  270.             static float gamma; /* peut être <0 */
  271.             static float c;     /* [0,PI]      */
  272.             static float a;     /* [0,PI]      */
  273.             static float aux1;
  274.             static float aux2;
  275.             /**********************************************************/
  276.             /* changement de repère en coordonnées sphériques         */
  277.             /* principe du triangle sphérique ABC: A est le pôle nord */
  278.             /* B parcours la carte, C est le nouveau pôle nord        */
  279.             /* on connait le coté b (soit [A,C]); il ne bouge pas, et */
  280.             /* calculé dans Update_New_Clip                           */
  281.             /* c dépend de la latitude du point B considéré           */
  282.             /* l'angle alpha (angle au sommet A) dépend de la longitu-*/
  283.             /* de du point B considéré                                */
  284.             /* Les formules de résolution du triangle sphérique nous  */
  285.             /* donnent le côté a et l'angle gamma, donnant respective-*/
  286.             /* ment les nouvelles latitude et longitude du point B    */
  287.             /**********************************************************/
  288.             alpha=(x-mx)*PI/18000;
  289.             c=(y+9000)*PI/18000;
  290.             if (Pref.prof==5 && FastsincosBase!=NULL)
  291.             {
  292.                 /* calcul de a et gamma                                   */
  293.                 aux1=cos_b2*fastcos(c)+sin_b2*fastsin(c)*fastcos(alpha);
  294.                 a=acos(aux1);
  295.                 aux2=fastsin(c)*fastsin(alpha)/fastsin(a);
  296.                 if (a!=0)
  297.                     /* asin nous donne la solution entre -PI/2 et PI/2 */
  298.                     gamma=asin(aux2);
  299.                 else
  300.                     gamma=0;
  301.  
  302.                 /* détermination des x et y*/
  303.                 y=-aux1*clip_max.maxy;
  304. /*                if (c>a && gamma<alpha)
  305.                     x=fastsin(a)*clip_max.maxx;
  306.                 else
  307.                 if (c<a && gamma>alpha)
  308.                     x=-sin(a)*clip_max.maxx;
  309.                 else*/
  310.                     /* la solution donnée par asin était bonne    */
  311. //                    if (gamma!=0 || x==mx || abs(x-mx)==18000)
  312.                         x=fastsin(a)*aux2*clip_max.maxx;
  313. /*                    else
  314.                         if ((x-mx+36000)%36000<18000)
  315.                             x=fastsin(a)*clip_max.maxx;
  316.                         else
  317.                             x=-fastsin(a)*clip_max.maxx;
  318. */            }
  319.             else
  320.             {
  321.                 /* calcul de a et gamma                                   */
  322.                 aux1=cos_b2*cos(c)+sin_b2*sin(c)*cos(alpha);
  323.                 a=acos(aux1);
  324.                 aux2=sin(c)*sin(alpha)/sin(a);
  325.                 if (a!=0)
  326.                     /* asin nous donne la solution entre -PI/2 et PI/2 */
  327.                     gamma=asin(aux2);
  328.                 else
  329.                     gamma=0;
  330.  
  331.                 /* détermination des x et y*/
  332.                 y=-aux1*clip_max.maxy;
  333.     /*          if (ABS(c)>ABS(a) && ABS(gamma)<ABS(alpha))
  334.                     x=sin(a)*clip_max.maxx;
  335.                 else
  336.                 if (ABS(c)<ABS(a) && ABS(gamma)>ABS(alpha))
  337.                     x=-sin(a)*clip_max.maxx;
  338.                 else */
  339.                     /* la solution donnée par asin était bonne    */
  340.                     if (gamma!=0 || x==mx || abs(x-mx)==18000)
  341.                         x=sin(a)*aux2*clip_max.maxx;
  342.                     else
  343.                         if ((x-mx+36000)%36000<18000)
  344.                             x=sin(a)*clip_max.maxx;
  345.                         else
  346.                             x=-sin(a)*clip_max.maxx;
  347.             }
  348.             /* transformation finale en coordonnées écran              */
  349. *px=(((x-Pref.clip_cur.minx+mx)*(Pref.win_height*win_ratio))/(Pref.clip_cur.maxx - Pref.clip_cur.minx));
  350. *py=(((y-Pref.clip_cur.miny+my)*(Pref.win_height))/(Pref.clip_cur.maxy - Pref.clip_cur.miny));
  351.         }
  352.             break;
  353.     case PROJ_FLAT:
  354.     default:
  355.         x=(((x-Pref.clip_cur.minx)*(Pref.win_width))/(Pref.clip_cur.maxx - Pref.clip_cur.minx));        
  356.         y=(((y-Pref.clip_cur.miny)*(Pref.win_height))/(Pref.clip_cur.maxy - Pref.clip_cur.miny));
  357.         *py=y;
  358.         *px=x;
  359.         break;
  360.     }
  361.     if (Pref.Flg_Proj_3D)
  362.     {
  363.             T3D_Convert(x,y,px,py);
  364.     }
  365. }
  366.  
  367. void 
  368. correct_clip(CLIP * clip)
  369. {
  370.         /*rend le clip rectangle*/
  371.         int sx=clip->maxx-clip->minx;
  372.         static int SX=0;
  373.         static int SY=0;
  374.         if (SX==0)
  375.         {
  376.                 SX=clip_max.maxx-clip_max.minx;
  377.                 SY=clip_max.maxy-clip_max.miny;
  378.         }
  379.         switch (Pref.Type_Proj)
  380.         {
  381.         case PROJ_GLOBE:
  382.                 clip->maxy=sx*SY/SX+clip->miny;
  383.                 break;  
  384.                 
  385.         case PROJ_MERCATOR:
  386.                 {
  387. /*              float coef=Map_Convert_Angle(SY*sx/SX)
  388.                                         +asin((double)(clip->miny/clip_max.maxy));*/
  389.                 double coef2=((double)(SY*sx))/(double)SX;
  390.                 double coef3=(double)coef2/(double)clip_max.maxy;
  391.                 double coef=coef3+sin(clip->miny*PI/18000);
  392.                 if (coef<=1)
  393.                         clip->maxy=asin(coef)*clip_max.maxy*2/PI;
  394.                 else
  395.                 {
  396.                         /*deuxieme cas*/
  397. /*
  398.                         clip->miny=sin(PI/2-Map_Convert_Angle(SY*sx/SX))
  399.                                                 *clip_max.maxy;
  400.                         clip->maxy=clip_max.maxy;*/
  401.                 }
  402.                 }
  403.                 break;
  404.         case PROJ_FLAT:
  405.                 clip->maxy=sx*SY/SX+clip->miny;
  406.                 break;
  407.         default:
  408.                 clip->maxy=sx*SY/SX+clip->miny;
  409.                 break;
  410.         }
  411. }                       
  412.         
  413.  
  414.