[1811] | 1 | //#define TEST_MOLLWEIDE
|
---|
| 2 | #include <stdlib.h>
|
---|
| 3 | #include <stdio.h>
|
---|
| 4 | #include <math.h>
|
---|
[2615] | 5 | #include "sopnamsp.h"
|
---|
[1811] | 6 | #include "mollweide.h"
|
---|
| 7 |
|
---|
| 8 | #if defined(TEST_MOLLWEIDE)
|
---|
| 9 | #include "timing.h"
|
---|
| 10 | double _t;
|
---|
| 11 | #endif
|
---|
| 12 |
|
---|
| 13 | /*! \ingroup XAstroPack
|
---|
| 14 | \brief La vraie projection Mollweide (longitude,latitude) -> (xmoll,ymoll)
|
---|
| 15 | \param input: longitude [0,2PI[ en radians
|
---|
| 16 | \param input: latitude [-Pi/2,Pi/2] en radians
|
---|
| 17 | \param output: xmoll dans [-2*sqrt(2),2*sqrt(2)]
|
---|
| 18 | \param output: ymoll dans [-sqrt(2),sqrt(2)]
|
---|
| 19 | \warning Probleme de precision pour t pres de Pi/2 (latitude pres de -Pi/2 ou Pi/2)
|
---|
[1813] | 20 | \warning Longitude=Pi=180 au milieu en xmoll=0
|
---|
[1811] | 21 | \return t
|
---|
| 22 | \verbatim
|
---|
| 23 | - La vraie MollWeide:
|
---|
| 24 | Equation implicite en t: 2*t+sin(2*t) = Pi*sin(lat) = a
|
---|
| 25 | | xmoll = 2*sqrt(2)*(lon-Pi)/Pi*cos(t)
|
---|
| 26 | | ymoll = sqrt(2)*sin(t)
|
---|
| 27 | --------------------------------------------------------
|
---|
| 28 | INPUT: longitude [0,2PI[ en radians
|
---|
| 29 | latitude [-Pi/2,Pi/2] en radians
|
---|
| 30 | OUTPUT: xmoll dans [-2*sqrt(2),2*sqrt(2)]
|
---|
| 31 | ymoll dans [-sqrt(2),sqrt(2)]
|
---|
| 32 | RETURN: t
|
---|
| 33 | --------------------------------------------------------
|
---|
| 34 | REMARQUES:
|
---|
| 35 | - xmoll varie [-2*sqrt(2),2*sqrt(2)]
|
---|
[1813] | 36 | centre galactique (longitude=0) au bord,
|
---|
| 37 | anti-centre (longitude=180) au milieu.
|
---|
| 38 | pour centre galactique (longitude=0) au milieu,
|
---|
| 39 | anti-centre (longitude=180) au bord:
|
---|
| 40 | xmoll' = 2*sqrt(2)*longitude/Pi*cos(t) = xmoll+2*sqrt(2)
|
---|
| 41 | xmoll' varie [0,4*sqrt(2)]
|
---|
[1811] | 42 | - ymoll varie [-sqrt(2),sqrt(2)] plan gal en 0
|
---|
| 43 | - t varie [-Pi/2,Pi/2]
|
---|
| 44 | - a varie [-Pi,Pi] (a <-> -a <==> t <-> -t)
|
---|
| 45 | - Pour latitude --> -latitude alors t --> -t
|
---|
| 46 | latitude=90 --> t=Pi/2
|
---|
| 47 | latitude=-90 --> t=-Pi/2
|
---|
| 48 | RAPIDITE: 10^6 appels a MollWeide_XY en 5.2 seconde (PC Linux 500Mhz)
|
---|
| 49 | PRECISION: eps=1e-16 -> |t-tvrai|<1e-6 toujours
|
---|
| 50 | |t-tvrai|<1e-10 pour t<1.57
|
---|
| 51 | |t-tvrai|<1e-10 pour t<1.57
|
---|
| 52 | |t-tvrai|<3e-14 pour t<1.50
|
---|
| 53 | |t-tvrai|<5e-15 pour t<1.40
|
---|
| 54 | \endverbatim
|
---|
| 55 | */
|
---|
| 56 | double MollWeide_XY(double longitude,double latitude,double* xmoll,double*ymoll)
|
---|
| 57 | {
|
---|
| 58 | int i,niter=2;
|
---|
| 59 | double tini,t0,f0,f1,f2,e1,e2,dt,a,asgn,eps=1.e-16;
|
---|
| 60 |
|
---|
| 61 | a = M_PI*sin(latitude);
|
---|
| 62 | if(a<0.) {a*=-1.; asgn=-1.;} else asgn=1.;
|
---|
| 63 | // premier guess t1 ok pour t->0 (droite)
|
---|
| 64 | f1 = a/4.;
|
---|
| 65 | // second guess t2 ok pour t->Pi/2 (DL en Pi/2, type a+b*x^3)
|
---|
| 66 | f2 = 0.75*(M_PI-a);
|
---|
| 67 | if(f2>0.) f2=M_PI_2-pow(f2,0.33333333333333333); else f2=M_PI_2;
|
---|
| 68 | e1 = 2.*f1+sin(2.*f1)-a; if(e1<0.) e1 *= -1.;
|
---|
| 69 | e2 = 2.*f2+sin(2.*f2)-a; if(e2<0.) e2 *= -1.;
|
---|
| 70 | tini = (f1*e2+f2*e1)/(e1+e2);
|
---|
| 71 | //printf("---> a=%e tini=%e f1=%e f2=%e\n",a,tini,f1,f2);
|
---|
| 72 | if(tini<0) tini=0.; else if(tini>M_PI_2) tini=M_PI_2;
|
---|
| 73 | // L'approximation par le DL au 2sd ordre, on resoud:
|
---|
| 74 | // a = f(x0) + f'(x0)*(x-x0) + f''(x0)*(x-x0)^2/2
|
---|
| 75 | t0 = tini;
|
---|
| 76 | for(i=0;i<niter;i++) {
|
---|
| 77 | e1 = sin(2.*t0);
|
---|
| 78 | f0 = 2.*t0+e1; // (f(x0)) >=0 sur 0<t<Pi/2
|
---|
| 79 | //printf(".%d a-f0=%e\n",i,a-f0);
|
---|
| 80 | if(a-f0<eps && a-f0>-eps) break;
|
---|
| 81 | f1 = 2.*(1.+cos(2.*t0)); // (f'(x0)) >=0 sur 0<t<Pi/2
|
---|
| 82 | f2 = -2.*e1; // (f''(x0)/2) <=0 sur 0<t<Pi/2
|
---|
| 83 | dt = f1*f1 - 4.*(f0-a)*f2;
|
---|
| 84 | //printf(".%d f0=%e f1=%e f2=%e disc=%e\n",i,f0,f1,f2,dt);
|
---|
| 85 | if(dt>=0.) dt=-f1+sqrt(dt); else dt=0.;
|
---|
| 86 | ////printf(".%d -b+sqrt(disc)=%e\n",i,dt);
|
---|
| 87 | if(f2<-eps && (dt>eps || dt<-eps)) { // approx parabolique
|
---|
| 88 | t0 += dt/(2.*f2);
|
---|
| 89 | //printf(".%d parab t0=%e\n",i,t0);
|
---|
| 90 | } else if(f1>eps) {
|
---|
| 91 | t0 += (a-f0)/f1; // approx lineaire
|
---|
| 92 | //printf(".%d lin t0=%e\n",i,t0);
|
---|
| 93 | }
|
---|
| 94 | if(t0<0) t0=0.; else if(t0>M_PI_2) t0=M_PI_2;
|
---|
| 95 | }
|
---|
| 96 |
|
---|
| 97 | t0 *= asgn;
|
---|
| 98 | *xmoll = 2.*M_SQRT2*(longitude-M_PI)/M_PI*cos(t0);
|
---|
| 99 | *ymoll = M_SQRT2*sin(t0);
|
---|
| 100 |
|
---|
| 101 | #if defined(TEST_MOLLWEIDE)
|
---|
| 102 | _t = tini;
|
---|
| 103 | #endif
|
---|
| 104 |
|
---|
| 105 | return t0;
|
---|
| 106 | }
|
---|
| 107 |
|
---|
| 108 | /*! \ingroup XAstroPack
|
---|
| 109 | \brief La vraie projection Mollweide (xmoll,ymoll) -> (longitude,latitude)
|
---|
| 110 | \param input: xmoll dans [-2*sqrt(2),2*sqrt(2)]
|
---|
| 111 | \param input: ymoll dans [-sqrt(2),sqrt(2)]
|
---|
| 112 | \param output: longitude [0,2PI[ en radians
|
---|
| 113 | \param output: latitude [-Pi/2,Pi/2] en radians
|
---|
| 114 | \return t
|
---|
| 115 | \warning Probleme de precision pour t pres de Pi/2 (latitude pres de -Pi/2 ou Pi/2)
|
---|
[1813] | 116 | \warning Longitude=Pi=180 au milieu en xmoll=0
|
---|
[1811] | 117 | \verbatim
|
---|
| 118 | --------------------------------------------------------
|
---|
| 119 | INPUT: xmoll dans [-2*sqrt(2),2*sqrt(2)]
|
---|
| 120 | ymoll dans [-sqrt(2),sqrt(2)]
|
---|
| 121 | OUTPUT: longitude [0,2PI[ en radians
|
---|
| 122 | latitude [-Pi/2,Pi/2] en radians
|
---|
| 123 | RETURN: t if OK
|
---|
| 124 | -991 : xmoll hors limites
|
---|
| 125 | -992 : ymoll hors limites
|
---|
| 126 | -993 : non inversible (poles de la sphere)
|
---|
| 127 | --------------------------------------------------------
|
---|
| 128 | \endverbatim
|
---|
| 129 | */
|
---|
| 130 | double MollWeide_LL(double xmoll,double ymoll,double* longitude,double* latitude)
|
---|
| 131 | {
|
---|
| 132 | double t,a,eps=1.e-20;
|
---|
| 133 |
|
---|
| 134 | *longitude = *latitude = -999.;
|
---|
| 135 | if(xmoll<-2*M_SQRT2 || xmoll>2*M_SQRT2) return(-901.);
|
---|
| 136 | if(ymoll<-M_SQRT2 || ymoll>M_SQRT2) return(-902.);
|
---|
| 137 |
|
---|
| 138 | t = ymoll/M_SQRT2; if(t<-1.) t=-1.; else if(t>1.) t=1.;
|
---|
| 139 | t = asin(t);
|
---|
| 140 |
|
---|
| 141 | a = (2.*t+sin(2.*t))/M_PI; if(a<-1.) a=-1.; else if(a>1.) a=1.;
|
---|
| 142 | *latitude = asin(a);
|
---|
| 143 |
|
---|
| 144 | a= cos(t);
|
---|
| 145 | if(fabs(a)<eps) return(-903.);
|
---|
| 146 | *longitude = M_PI*(xmoll/(2.*M_SQRT2*a)+1.);
|
---|
| 147 |
|
---|
| 148 | return(t);
|
---|
| 149 | }
|
---|
| 150 |
|
---|
| 151 | ///////////////////////////////////////
|
---|
| 152 | ////////// Programme de test //////////
|
---|
| 153 | ///////////////////////////////////////
|
---|
| 154 | #if defined(TEST_MOLLWEIDE)
|
---|
| 155 | int main()
|
---|
| 156 | {
|
---|
| 157 | FILE *file = fopen("mollw.data","w");
|
---|
| 158 | FILE *file2 = fopen("mollw2.data","w");
|
---|
| 159 | int i,npoints=250000;
|
---|
| 160 | double tvrai,lon,lat,x,y,tcalc,a;
|
---|
| 161 | double longitude,latitude,tll,txy;
|
---|
| 162 | InitTim();
|
---|
| 163 |
|
---|
| 164 | tll=txy=0.;
|
---|
| 165 | for(i=0;i<=npoints;i++) {
|
---|
| 166 | //tvrai = 1.57078;
|
---|
| 167 | tvrai = (double)i*(M_PI_2) /(double)npoints;
|
---|
| 168 | a = 2.*tvrai+sin(2.*tvrai);
|
---|
| 169 | lat = asin(a/M_PI);
|
---|
| 170 | tcalc = MollWeide_XY(0.,lat,&x,&y);
|
---|
| 171 | if(fabs(tcalc-tvrai)>tll) {tll = fabs(tcalc-tvrai); txy=tvrai;}
|
---|
| 172 | if(file) fprintf(file,"%e %e %e %e %e\n",tvrai,tcalc,_t,a,tcalc-tvrai);
|
---|
| 173 | }
|
---|
| 174 | printf("Maximum difference: %e for t=%.10e=Pi/2-%.10e\n",tll,txy,M_PI_2-txy);
|
---|
| 175 | if(file) fclose(file);
|
---|
| 176 | PrtTim("--- End of first testing ---");
|
---|
| 177 |
|
---|
| 178 | for(lon=0.;lon<360.;lon+=2.5) for(lat=-90.;lat<=90.1;lat+=2.5) {
|
---|
| 179 | txy = MollWeide_XY(lon*M_PI/180.,lat*M_PI/180.,&x,&y);
|
---|
| 180 | tll = MollWeide_LL(x,y,&longitude,&latitude);
|
---|
| 181 | longitude *= 180./M_PI;
|
---|
| 182 | latitude *= 180./M_PI;
|
---|
| 183 | if(file2) fprintf(file2,"%e %e %e %e %e %e %e %e %e %e\n"
|
---|
| 184 | ,lon,lat,x,y,longitude,latitude,txy,tll,longitude-lon,latitude-lat);
|
---|
| 185 | }
|
---|
| 186 | if(file2) fclose(file2);
|
---|
| 187 | PrtTim("--- End of second testing ---");
|
---|
| 188 | return 0;
|
---|
| 189 | }
|
---|
| 190 | #endif
|
---|