| [2551] | 1 | /* DoD NIMA World Magnetic Model.
|
|---|
| 2 | * from http://www.ngdc.noaa.gov
|
|---|
| 3 | *
|
|---|
| 4 | #define TEST_MAIN
|
|---|
| 5 | */
|
|---|
| 6 |
|
|---|
| 7 |
|
|---|
| 8 | #include <math.h>
|
|---|
| 9 | #include <stdio.h>
|
|---|
| 10 | #include <string.h>
|
|---|
| 11 | #include <errno.h>
|
|---|
| 12 |
|
|---|
| 13 | #include "astro.h"
|
|---|
| 14 |
|
|---|
| 15 | static char mfn[] = "wmm.cof"; /* file with model coefficients */
|
|---|
| 16 |
|
|---|
| 17 | static int geomag(FILE *wmmdat, int *maxdeg);
|
|---|
| 18 | static int geomg1(FILE *wmmdat, float alt, float glat, float glon,
|
|---|
| 19 | float t, float *dec, float *mdp, float *ti, float *gv);
|
|---|
| 20 |
|
|---|
| 21 | /* compute magnetic declination for given location, elevation and time.
|
|---|
| 22 | * sign is such that mag bearing = true az + mag deviation.
|
|---|
| 23 | * return 0 if ok, -1 if no model file, -2 if time outside model range.
|
|---|
| 24 | * fill err[] with excuse if return < 0.
|
|---|
| 25 | */
|
|---|
| 26 | int
|
|---|
| 27 | magdecl (
|
|---|
| 28 | double l, double L, /* geodesic lat, +N, long, +E, rads */
|
|---|
| 29 | double e, /* elevation, m */
|
|---|
| 30 | double y, /* time, decimal year */
|
|---|
| 31 | char *dir, /* dir for model file */
|
|---|
| 32 | double *mdp, /* magnetic deviation, rads E of N */
|
|---|
| 33 | char *err) /* err message if return < 0 */
|
|---|
| 34 | {
|
|---|
| 35 | float dlat = raddeg(l);
|
|---|
| 36 | float dlon = raddeg(L);
|
|---|
| 37 | float alt = e/1000.;
|
|---|
| 38 | int maxdeg = 12;
|
|---|
| 39 | float dec, dp, ti, gv;
|
|---|
| 40 | char mfile[1024];
|
|---|
| 41 | FILE *wmmdat;
|
|---|
| 42 | int s;
|
|---|
| 43 |
|
|---|
| 44 | /* open model file */
|
|---|
| 45 | sprintf (mfile, "%s/%s", dir, mfn);
|
|---|
| 46 | wmmdat = fopen (mfile, "r");
|
|---|
| 47 | if (!wmmdat) {
|
|---|
| 48 | sprintf (err, "%s: %s", mfile, strerror(errno));
|
|---|
| 49 | return (-1);
|
|---|
| 50 | }
|
|---|
| 51 |
|
|---|
| 52 | /* compute deviation */
|
|---|
| 53 | geomag(wmmdat, &maxdeg);
|
|---|
| 54 | s = geomg1(wmmdat,alt,dlat,dlon,y,&dec,&dp,&ti,&gv);
|
|---|
| 55 | fclose(wmmdat);
|
|---|
| 56 | if (s < 0) {
|
|---|
| [2641] | 57 | sprintf (err, "%s: Magnetic model only available for %g .. %g. See http://www.ngdc.noaa.gov", mfile, ti, ti+5);
|
|---|
| [2551] | 58 | return (-2);
|
|---|
| 59 | }
|
|---|
| 60 | *mdp = degrad(dec);
|
|---|
| 61 | return (0);
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 | #if defined(TEST_MAIN)
|
|---|
| 65 |
|
|---|
| 66 | int
|
|---|
| 67 | main(int ac, char *av[])
|
|---|
| 68 | {
|
|---|
| 69 | char err[1024];
|
|---|
| 70 | float altm, dlat, dlon;
|
|---|
| 71 | float t;
|
|---|
| 72 | double dec;
|
|---|
| 73 |
|
|---|
| 74 | S1:
|
|---|
| 75 | printf("\n\n\n ENTER LATITUDE IN DECIMAL DEGREES (+25.0)\n");
|
|---|
| 76 | scanf("%f", &dlat);
|
|---|
| 77 |
|
|---|
| 78 | printf(" ENTER LONGITUDE IN DECIMAL DEGREES (-100.0)\n");
|
|---|
| 79 | scanf("%f", &dlon);
|
|---|
| 80 |
|
|---|
| 81 | printf(" ENTER ALTITUDE IN METERS\n");
|
|---|
| 82 | scanf("%f", &altm);
|
|---|
| 83 |
|
|---|
| 84 | printf(" ENTER TIME IN DECIMAL YEAR\n");
|
|---|
| 85 | scanf("%f",&t);
|
|---|
| 86 |
|
|---|
| 87 | if (magdecl (degrad(dlat), degrad(dlon), altm, t, "auxil", &dec,
|
|---|
| 88 | err) < 0) {
|
|---|
| 89 | printf ("%s\n", err);
|
|---|
| 90 | return(1);
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | printf("\n LATITUDE: = %-7.2f DEG",dlat);
|
|---|
| 94 | printf("\n LONGITUDE: = %-7.2f DEG\n",dlon);
|
|---|
| 95 | printf("\n ALTITUDE = %.2f METERS",altm);
|
|---|
| 96 | printf("\n DATE = %-5.1f\n",t);
|
|---|
| 97 |
|
|---|
| 98 | printf("\n\t\t\t OUTPUT\n\t\t\t ------");
|
|---|
| 99 |
|
|---|
| 100 | printf("\n DEC = %-7.2f DEG", raddeg(dec));
|
|---|
| 101 |
|
|---|
| 102 | printf("\n\n\n DO YOU NEED MORE POINT DATA? (y or n)\n");
|
|---|
| 103 | scanf("%s", err);
|
|---|
| 104 | if ((err[0] =='y')||(err[0] == 'Y')) goto S1;
|
|---|
| 105 |
|
|---|
| 106 | return(0);
|
|---|
| 107 | }
|
|---|
| 108 | #endif /* defined(TEST_MAIN) */
|
|---|
| 109 |
|
|---|
| 110 | /*************************************************************************
|
|---|
| 111 | * return 0 if ok, -1 if time is out of range with base epoc in *ti
|
|---|
| 112 | */
|
|---|
| 113 |
|
|---|
| 114 | static int E0000(FILE *wmmdat, int IENTRY, int *maxdeg, float alt,
|
|---|
| 115 | float glat, float glon, float t, float *dec, float *mdp, float *ti,
|
|---|
| 116 | float *gv)
|
|---|
| 117 | {
|
|---|
| 118 | static int maxord,i,icomp,n,m,j,D1,D2,D3,D4;
|
|---|
| 119 | static float c[13][13],cd[13][13],tc[13][13],dp[13][13],snorm[169],
|
|---|
| 120 | sp[13],cp[13],fn[13],fm[13],pp[13],k[13][13],pi,dtr,a,b,re,
|
|---|
| 121 | a2,b2,c2,a4,b4,c4,epoc,gnm,hnm,dgnm,dhnm,flnmj,otime,oalt,
|
|---|
| 122 | olat,olon,dt,rlon,rlat,srlon,srlat,crlon,crlat,srlat2,
|
|---|
| 123 | crlat2,q,q1,q2,ct,st,r2,r,d,ca,sa,aor,ar,br,bt,bp,bpp,
|
|---|
| 124 | par,temp1,temp2,parp,bx,by,bz,bh;
|
|---|
| 125 | static char model[20], c_str[81], c_new[5];
|
|---|
| 126 | static float *p = snorm;
|
|---|
| 127 |
|
|---|
| 128 | switch(IENTRY){case 0: goto GEOMAG; case 1: goto GEOMG1;}
|
|---|
| 129 |
|
|---|
| 130 | GEOMAG:
|
|---|
| 131 |
|
|---|
| 132 | /* INITIALIZE CONSTANTS */
|
|---|
| 133 | maxord = *maxdeg;
|
|---|
| 134 | sp[0] = 0.0;
|
|---|
| 135 | cp[0] = *p = pp[0] = 1.0;
|
|---|
| 136 | dp[0][0] = 0.0;
|
|---|
| 137 | a = 6378.137;
|
|---|
| 138 | b = 6356.7523142;
|
|---|
| 139 | re = 6371.2;
|
|---|
| 140 | a2 = a*a;
|
|---|
| 141 | b2 = b*b;
|
|---|
| 142 | c2 = a2-b2;
|
|---|
| 143 | a4 = a2*a2;
|
|---|
| 144 | b4 = b2*b2;
|
|---|
| 145 | c4 = a4 - b4;
|
|---|
| 146 |
|
|---|
| 147 | /* READ WORLD MAGNETIC MODEL SPHERICAL HARMONIC COEFFICIENTS */
|
|---|
| 148 | c[0][0] = 0.0;
|
|---|
| 149 | cd[0][0] = 0.0;
|
|---|
| 150 | fgets(c_str, 80, wmmdat);
|
|---|
| 151 | sscanf(c_str,"%f%s",&epoc,model);
|
|---|
| 152 | S3:
|
|---|
| 153 | fgets(c_str, 80, wmmdat);
|
|---|
| 154 | /* CHECK FOR LAST LINE IN FILE */
|
|---|
| 155 | for (i=0; i<4 && (c_str[i] != '\0'); i++)
|
|---|
| 156 | {
|
|---|
| 157 | c_new[i] = c_str[i];
|
|---|
| 158 | c_new[i+1] = '\0';
|
|---|
| 159 | }
|
|---|
| 160 | icomp = strcmp("9999", c_new);
|
|---|
| 161 | if (icomp == 0) goto S4;
|
|---|
| 162 | /* END OF FILE NOT ENCOUNTERED, GET VALUES */
|
|---|
| 163 | sscanf(c_str,"%d%d%f%f%f%f",&n,&m,&gnm,&hnm,&dgnm,&dhnm);
|
|---|
| 164 | if (m <= n)
|
|---|
| 165 | {
|
|---|
| 166 | c[m][n] = gnm;
|
|---|
| 167 | cd[m][n] = dgnm;
|
|---|
| 168 | if (m != 0)
|
|---|
| 169 | {
|
|---|
| 170 | c[n][m-1] = hnm;
|
|---|
| 171 | cd[n][m-1] = dhnm;
|
|---|
| 172 | }
|
|---|
| 173 | }
|
|---|
| 174 | goto S3;
|
|---|
| 175 |
|
|---|
| 176 | /* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */
|
|---|
| 177 | S4:
|
|---|
| 178 | *snorm = 1.0;
|
|---|
| 179 | for (n=1; n<=maxord; n++)
|
|---|
| 180 | {
|
|---|
| 181 | *(snorm+n) = *(snorm+n-1)*(float)(2*n-1)/(float)n;
|
|---|
| 182 | j = 2;
|
|---|
| 183 | for (m=0,D1=1,D2=(n-m+D1)/D1; D2>0; D2--,m+=D1)
|
|---|
| 184 | {
|
|---|
| 185 | k[m][n] = (float)(((n-1)*(n-1))-(m*m))/(float)((2*n-1)*(2*n-3));
|
|---|
| 186 | if (m > 0)
|
|---|
| 187 | {
|
|---|
| 188 | flnmj = (float)((n-m+1)*j)/(float)(n+m);
|
|---|
| 189 | *(snorm+n+m*13) = *(snorm+n+(m-1)*13)*sqrt(flnmj);
|
|---|
| 190 | j = 1;
|
|---|
| 191 | c[n][m-1] = *(snorm+n+m*13)*c[n][m-1];
|
|---|
| 192 | cd[n][m-1] = *(snorm+n+m*13)*cd[n][m-1];
|
|---|
| 193 | }
|
|---|
| 194 | c[m][n] = *(snorm+n+m*13)*c[m][n];
|
|---|
| 195 | cd[m][n] = *(snorm+n+m*13)*cd[m][n];
|
|---|
| 196 | }
|
|---|
| 197 | fn[n] = (float)(n+1);
|
|---|
| 198 | fm[n] = (float)n;
|
|---|
| 199 | }
|
|---|
| 200 | k[1][1] = 0.0;
|
|---|
| 201 |
|
|---|
| 202 | otime = oalt = olat = olon = -1000.0;
|
|---|
| 203 | return (0);
|
|---|
| 204 |
|
|---|
| 205 | /*************************************************************************/
|
|---|
| 206 |
|
|---|
| 207 | GEOMG1:
|
|---|
| 208 |
|
|---|
| 209 | dt = t - epoc;
|
|---|
| 210 | if (otime < 0.0 && (dt < 0.0 || dt > 5.0)) {
|
|---|
| 211 | *ti = epoc; /* pass back base time for diag msg */
|
|---|
| 212 | return (-1);
|
|---|
| 213 | }
|
|---|
| 214 |
|
|---|
| 215 | pi = 3.14159265359;
|
|---|
| 216 | dtr = pi/180.0;
|
|---|
| 217 | rlon = glon*dtr;
|
|---|
| 218 | rlat = glat*dtr;
|
|---|
| 219 | srlon = sin(rlon);
|
|---|
| 220 | srlat = sin(rlat);
|
|---|
| 221 | crlon = cos(rlon);
|
|---|
| 222 | crlat = cos(rlat);
|
|---|
| 223 | srlat2 = srlat*srlat;
|
|---|
| 224 | crlat2 = crlat*crlat;
|
|---|
| 225 | sp[1] = srlon;
|
|---|
| 226 | cp[1] = crlon;
|
|---|
| 227 |
|
|---|
| 228 | /* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */
|
|---|
| 229 | if (alt != oalt || glat != olat)
|
|---|
| 230 | {
|
|---|
| 231 | q = sqrt(a2-c2*srlat2);
|
|---|
| 232 | q1 = alt*q;
|
|---|
| 233 | q2 = ((q1+a2)/(q1+b2))*((q1+a2)/(q1+b2));
|
|---|
| 234 | ct = srlat/sqrt(q2*crlat2+srlat2);
|
|---|
| 235 | st = sqrt(1.0-(ct*ct));
|
|---|
| 236 | r2 = (alt*alt)+2.0*q1+(a4-c4*srlat2)/(q*q);
|
|---|
| 237 | r = sqrt(r2);
|
|---|
| 238 | d = sqrt(a2*crlat2+b2*srlat2);
|
|---|
| 239 | ca = (alt+d)/r;
|
|---|
| 240 | sa = c2*crlat*srlat/(r*d);
|
|---|
| 241 | }
|
|---|
| 242 | if (glon != olon)
|
|---|
| 243 | {
|
|---|
| 244 | for (m=2; m<=maxord; m++)
|
|---|
| 245 | {
|
|---|
| 246 | sp[m] = sp[1]*cp[m-1]+cp[1]*sp[m-1];
|
|---|
| 247 | cp[m] = cp[1]*cp[m-1]-sp[1]*sp[m-1];
|
|---|
| 248 | }
|
|---|
| 249 | }
|
|---|
| 250 | aor = re/r;
|
|---|
| 251 | ar = aor*aor;
|
|---|
| 252 | br = bt = bp = bpp = 0.0;
|
|---|
| 253 | for (n=1; n<=maxord; n++)
|
|---|
| 254 | {
|
|---|
| 255 | ar = ar*aor;
|
|---|
| 256 | for (m=0,D3=1,D4=(n+m+D3)/D3; D4>0; D4--,m+=D3)
|
|---|
| 257 | {
|
|---|
| 258 | /*
|
|---|
| 259 | COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
|
|---|
| 260 | AND DERIVATIVES VIA RECURSION RELATIONS
|
|---|
| 261 | */
|
|---|
| 262 | if (alt != oalt || glat != olat)
|
|---|
| 263 | {
|
|---|
| 264 | if (n == m)
|
|---|
| 265 | {
|
|---|
| 266 | *(p+n+m*13) = st**(p+n-1+(m-1)*13);
|
|---|
| 267 | dp[m][n] = st*dp[m-1][n-1]+ct**(p+n-1+(m-1)*13);
|
|---|
| 268 | goto S50;
|
|---|
| 269 | }
|
|---|
| 270 | if (n == 1 && m == 0)
|
|---|
| 271 | {
|
|---|
| 272 | *(p+n+m*13) = ct**(p+n-1+m*13);
|
|---|
| 273 | dp[m][n] = ct*dp[m][n-1]-st**(p+n-1+m*13);
|
|---|
| 274 | goto S50;
|
|---|
| 275 | }
|
|---|
| 276 | if (n > 1 && n != m)
|
|---|
| 277 | {
|
|---|
| 278 | if (m > n-2) *(p+n-2+m*13) = 0.0;
|
|---|
| 279 | if (m > n-2) dp[m][n-2] = 0.0;
|
|---|
| 280 | *(p+n+m*13) = ct**(p+n-1+m*13)-k[m][n]**(p+n-2+m*13);
|
|---|
| 281 | dp[m][n] = ct*dp[m][n-1] - st**(p+n-1+m*13)-k[m][n]*dp[m][n-2];
|
|---|
| 282 | }
|
|---|
| 283 | }
|
|---|
| 284 | S50:
|
|---|
| 285 | /*
|
|---|
| 286 | TIME ADJUST THE GAUSS COEFFICIENTS
|
|---|
| 287 | */
|
|---|
| 288 | if (t != otime)
|
|---|
| 289 | {
|
|---|
| 290 | tc[m][n] = c[m][n]+dt*cd[m][n];
|
|---|
| 291 | if (m != 0) tc[n][m-1] = c[n][m-1]+dt*cd[n][m-1];
|
|---|
| 292 | }
|
|---|
| 293 | /*
|
|---|
| 294 | ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
|
|---|
| 295 | */
|
|---|
| 296 | par = ar**(p+n+m*13);
|
|---|
| 297 | if (m == 0)
|
|---|
| 298 | {
|
|---|
| 299 | temp1 = tc[m][n]*cp[m];
|
|---|
| 300 | temp2 = tc[m][n]*sp[m];
|
|---|
| 301 | }
|
|---|
| 302 | else
|
|---|
| 303 | {
|
|---|
| 304 | temp1 = tc[m][n]*cp[m]+tc[n][m-1]*sp[m];
|
|---|
| 305 | temp2 = tc[m][n]*sp[m]-tc[n][m-1]*cp[m];
|
|---|
| 306 | }
|
|---|
| 307 | bt = bt-ar*temp1*dp[m][n];
|
|---|
| 308 | bp += (fm[m]*temp2*par);
|
|---|
| 309 | br += (fn[n]*temp1*par);
|
|---|
| 310 | /*
|
|---|
| 311 | SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES
|
|---|
| 312 | */
|
|---|
| 313 | if (st == 0.0 && m == 1)
|
|---|
| 314 | {
|
|---|
| 315 | if (n == 1) pp[n] = pp[n-1];
|
|---|
| 316 | else pp[n] = ct*pp[n-1]-k[m][n]*pp[n-2];
|
|---|
| 317 | parp = ar*pp[n];
|
|---|
| 318 | bpp += (fm[m]*temp2*parp);
|
|---|
| 319 | }
|
|---|
| 320 | }
|
|---|
| 321 | }
|
|---|
| 322 | if (st == 0.0) bp = bpp;
|
|---|
| 323 | else bp /= st;
|
|---|
| 324 | /*
|
|---|
| 325 | ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO
|
|---|
| 326 | GEODETIC COORDINATES
|
|---|
| 327 | */
|
|---|
| 328 | bx = -bt*ca-br*sa;
|
|---|
| 329 | by = bp;
|
|---|
| 330 | bz = bt*sa-br*ca;
|
|---|
| 331 | /*
|
|---|
| 332 | COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND
|
|---|
| 333 | TOTAL INTENSITY (TI)
|
|---|
| 334 | */
|
|---|
| 335 | bh = sqrt((bx*bx)+(by*by));
|
|---|
| 336 | *ti = sqrt((bh*bh)+(bz*bz));
|
|---|
| 337 | *dec = atan2(by,bx)/dtr;
|
|---|
| 338 | *mdp = atan2(bz,bh)/dtr;
|
|---|
| 339 | /*
|
|---|
| 340 | COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT
|
|---|
| 341 | GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC
|
|---|
| 342 | (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES)
|
|---|
| 343 |
|
|---|
| 344 | OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0
|
|---|
| 345 | */
|
|---|
| 346 | *gv = -999.0;
|
|---|
| 347 | if (fabs(glat) >= 55.)
|
|---|
| 348 | {
|
|---|
| 349 | if (glat > 0.0 && glon >= 0.0) *gv = *dec-glon;
|
|---|
| 350 | if (glat > 0.0 && glon < 0.0) *gv = *dec+fabs(glon);
|
|---|
| 351 | if (glat < 0.0 && glon >= 0.0) *gv = *dec+glon;
|
|---|
| 352 | if (glat < 0.0 && glon < 0.0) *gv = *dec-fabs(glon);
|
|---|
| 353 | if (*gv > +180.0) *gv -= 360.0;
|
|---|
| 354 | if (*gv < -180.0) *gv += 360.0;
|
|---|
| 355 | }
|
|---|
| 356 | otime = t;
|
|---|
| 357 | oalt = alt;
|
|---|
| 358 | olat = glat;
|
|---|
| 359 | olon = glon;
|
|---|
| 360 | return (0);
|
|---|
| 361 | }
|
|---|
| 362 |
|
|---|
| 363 | /*************************************************************************/
|
|---|
| 364 |
|
|---|
| 365 | static int
|
|---|
| 366 | geomag(FILE *wmmdat, int *maxdeg)
|
|---|
| 367 | {
|
|---|
| 368 | return (E0000(wmmdat,0,maxdeg,0.0,0.0,0.0,0.0,NULL,NULL,NULL,NULL));
|
|---|
| 369 | }
|
|---|
| 370 |
|
|---|
| 371 | /*************************************************************************/
|
|---|
| 372 |
|
|---|
| 373 | static int
|
|---|
| 374 | geomg1(FILE *wmmdat, float alt, float glat, float glon, float t,
|
|---|
| 375 | float *dec, float *mdp, float *ti, float *gv)
|
|---|
| 376 | {
|
|---|
| 377 | return (E0000(wmmdat,1,NULL,alt,glat,glon,t,dec,mdp,ti,gv));
|
|---|
| 378 | }
|
|---|
| 379 |
|
|---|
| 380 | /* For RCS Only -- Do Not Edit */
|
|---|
| [3654] | 381 | static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: magdecl.c,v $ $Date: 2009-07-16 10:34:37 $ $Revision: 1.7 $ $Name: not supported by cvs2svn $"};
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
