OpenCPN Partial API docs
Loading...
Searching...
No Matches
GeomagnetismLibrary.c
1
2#include <stdio.h>
3#include <string.h>
4#include <math.h>
5#include <stdlib.h>
6#include <ctype.h>
7#include <assert.h>
8#include "GeomagnetismHeader.h"
9
10#ifndef OPENCPN
11#define OPENCPN 1
12#endif
13
14/* $Id: GeomagnetismLibrary.c 1521 2017-01-24 17:52:41Z awoods $
15 *
16 * ABSTRACT
17 *
18 * The purpose of Geomagnetism Library is primarily to support the World
19 Magnetic Model (WMM) 2015-2020.
20 * It however is built to be used for spherical harmonic models of the Earth's
21 magnetic field
22 * generally and supports models even with a large (>>12) number of degrees. It
23 is also used in many
24 * other geomagnetic models distributed by NCEI.
25 *
26 * REUSE NOTES
27 *
28 * Geomagnetism Library is intended for reuse by any application that requires
29 * Computation of Geomagnetic field from a spherical harmonic model.
30 *
31 * REFERENCES
32 *
33 * Further information on Geoid can be found in the WMM Technical Documents.
34 *
35 *
36 * LICENSES
37 *
38 * The WMM source code is in the public domain and not licensed or under
39 copyright.
40 * The information and software may be used freely by the public. As required
41 by 17 U.S.C. 403,
42 * third parties producing copyrighted works consisting predominantly of the
43 material produced by
44 * U.S. government agencies must provide notice with such work(s) identifying
45 the U.S. Government material
46 * incorporated and stating that such material is not subject to copyright
47 protection.
48 *
49 * RESTRICTIONS
50 *
51 * Geomagnetism library has no restrictions.
52 *
53 * ENVIRONMENT
54 *
55 * Geomagnetism library was tested in the following environments
56 *
57 * 1. Red Hat Linux with GCC Compiler
58 * 2. MS Windows 7 with MinGW compiler
59 * 3. Sun Solaris with GCC Compiler
60 *
61 *
62
63
64 * National Centers for Environmental Information
65 * NOAA E/NE42, 325 Broadway
66 * Boulder, CO 80305 USA
67 * Attn: Arnaud Chulliat
68 * Phone: (303) 497-6522
69 * Email: Arnaud.Chulliat@noaa.gov
70
71 * Software and Model Support
72 * National Centers for Environmental Information
73 * NOAA E/NE42
74 * 325 Broadway
75 * Boulder, CO 80305 USA
76 * Attn: Adam Woods or Manoj Nair
77 * Phone: (303) 497-6640 or -4642
78 * Email: geomag.models@noaa.gov
79 * URL: http://www.ngdc.noaa.gov/Geomagnetic/WMM/DoDWMM.shtml
80
81
82 * For more details on the subroutines, please consult the WMM
83 * Technical Documentations at
84 * http://www.ngdc.noaa.gov/Geomagnetic/WMM/DoDWMM.shtml
85
86 * Nov 23, 2009
87 * Written by Manoj C Nair and Adam Woods
88 * Manoj.C.Nair@noaa.Gov
89 * Adam.Woods@noaa.gov
90 */
91
92/******************************************************************************
93 ************************************Wrapper***********************************
94 * This grouping consists of functions call groups of other functions to do a
95 * complete calculation of some sort. For example, the MAG_Geomag function
96 * does everything necessary to compute the geomagnetic elements from a given
97 * geodetic point in space and magnetic model adjusted for the appropriate
98 * date. These functions are the external functions necessary to create a
99 * program that uses or calculates the magnetic field.
100 ******************************************************************************
101 ******************************************************************************/
102
103int MAG_Geomag(MAGtype_Ellipsoid Ellip, MAGtype_CoordSpherical CoordSpherical,
104 MAGtype_CoordGeodetic CoordGeodetic,
105 MAGtype_MagneticModel *TimedMagneticModel,
106 MAGtype_GeoMagneticElements *GeoMagneticElements)
107/*
108The main subroutine that calls a sequence of WMM sub-functions to calculate the
109magnetic field elements for a single point. The function expects the model
110coefficients and point coordinates as input and returns the magnetic field
111elements and their rate of change. Though, this subroutine can be called
112successively to calculate a time series, profile or grid of magnetic field,
113these are better achieved by the subroutine MAG_Grid.
114
115INPUT: Ellip
116 CoordSpherical
117 CoordGeodetic
118 TimedMagneticModel
119
120OUTPUT : GeoMagneticElements
121
122CALLS: MAG_AllocateLegendreFunctionMemory(NumTerms); ( For storing the ALF
123functions ) MAG_ComputeSphericalHarmonicVariables( Ellip, CoordSpherical,
124TimedMagneticModel->nMax, &SphVariables); (Compute Spherical Harmonic variables
125) MAG_AssociatedLegendreFunction(CoordSpherical, TimedMagneticModel->nMax,
126LegendreFunction); Compute ALF MAG_Summation(LegendreFunction,
127TimedMagneticModel, SphVariables, CoordSpherical, &MagneticResultsSph);
128Accumulate the spherical harmonic coefficients
129 MAG_SecVarSummation(LegendreFunction, TimedMagneticModel,
130SphVariables, CoordSpherical, &MagneticResultsSphVar); Sum the Secular Variation
131Coefficients MAG_RotateMagneticVector(CoordSpherical, CoordGeodetic,
132MagneticResultsSph, &MagneticResultsGeo); Map the computed Magnetic fields to
133Geodetic coordinates MAG_CalculateGeoMagneticElements(&MagneticResultsGeo,
134GeoMagneticElements); Calculate the Geomagnetic elements
135 MAG_CalculateSecularVariationElements(MagneticResultsGeoVar,
136GeoMagneticElements); Calculate the secular variation of each of the Geomagnetic
137elements
138
139 */
140{
141 MAGtype_LegendreFunction *LegendreFunction;
142 MAGtype_SphericalHarmonicVariables *SphVariables;
143 int NumTerms;
144 MAGtype_MagneticResults MagneticResultsSph, MagneticResultsGeo,
145 MagneticResultsSphVar, MagneticResultsGeoVar;
146
147 NumTerms =
148 ((TimedMagneticModel->nMax + 1) * (TimedMagneticModel->nMax + 2) / 2);
149 LegendreFunction = MAG_AllocateLegendreFunctionMemory(
150 NumTerms); /* For storing the ALF functions */
151 SphVariables = MAG_AllocateSphVarMemory(TimedMagneticModel->nMax);
152 MAG_ComputeSphericalHarmonicVariables(
153 Ellip, CoordSpherical, TimedMagneticModel->nMax,
154 SphVariables); /* Compute Spherical Harmonic variables */
155 MAG_AssociatedLegendreFunction(CoordSpherical, TimedMagneticModel->nMax,
156 LegendreFunction); /* Compute ALF */
157 MAG_Summation(
158 LegendreFunction, TimedMagneticModel, *SphVariables, CoordSpherical,
159 &MagneticResultsSph); /* Accumulate the spherical harmonic coefficients*/
160 MAG_SecVarSummation(
161 LegendreFunction, TimedMagneticModel, *SphVariables, CoordSpherical,
162 &MagneticResultsSphVar); /*Sum the Secular Variation Coefficients */
163 MAG_RotateMagneticVector(
164 CoordSpherical, CoordGeodetic, MagneticResultsSph,
165 &MagneticResultsGeo); /* Map the computed Magnetic fields to Geodeitic
166 coordinates */
167 MAG_RotateMagneticVector(
168 CoordSpherical, CoordGeodetic, MagneticResultsSphVar,
169 &MagneticResultsGeoVar); /* Map the secular variation field components to
170 Geodetic coordinates*/
171 MAG_CalculateGeoMagneticElements(
172 &MagneticResultsGeo,
173 GeoMagneticElements); /* Calculate the Geomagnetic elements, Equation 19 ,
174 WMM Technical report */
175 MAG_CalculateSecularVariationElements(
176 MagneticResultsGeoVar,
177 GeoMagneticElements); /*Calculate the secular variation of each of the
178 Geomagnetic elements*/
179
180 MAG_FreeLegendreMemory(LegendreFunction);
181 MAG_FreeSphVarMemory(SphVariables);
182
183 return TRUE;
184} /*MAG_Geomag*/
185
186void MAG_Gradient(MAGtype_Ellipsoid Ellip, MAGtype_CoordGeodetic CoordGeodetic,
187 MAGtype_MagneticModel *TimedMagneticModel,
188 MAGtype_Gradient *Gradient) {
189 /*It should be noted that the x[2], y[2], and z[2] variables are NOT the same
190 coordinate system as the directions in which the gradients are taken. These
191 variables represent a Cartesian coordinate system where the Earth's center is
192 the origin, 'z' points up toward the North (rotational) pole and 'x' points
193 toward the prime meridian. 'y' points toward longitude = 90 degrees East.
194 The gradient is preformed along a local Cartesian coordinate system with the
195 origin at CoordGeodetic. 'z' points down toward the Earth's core, x points
196 North, tangent to the local longitude line, and 'y' points East, tangent to
197 the local latitude line.*/
198 double phiDelta = 0.01, /*DeltaY = 0.01,*/ hDelta = -1, x[2], y[2], z[2],
199 distance;
200
201 MAGtype_CoordSpherical AdjCoordSpherical;
202 MAGtype_CoordGeodetic AdjCoordGeodetic;
203 MAGtype_GeoMagneticElements GeomagneticElements, AdjGeoMagneticElements[2];
204
205 /*Initialization*/
206 MAG_GeodeticToSpherical(Ellip, CoordGeodetic, &AdjCoordSpherical);
207 MAG_Geomag(Ellip, AdjCoordSpherical, CoordGeodetic, TimedMagneticModel,
208 &GeomagneticElements);
209 AdjCoordGeodetic = MAG_CoordGeodeticAssign(CoordGeodetic);
210
211 /*Gradient along x*/
212
213 AdjCoordGeodetic.phi = CoordGeodetic.phi + phiDelta;
214 MAG_GeodeticToSpherical(Ellip, AdjCoordGeodetic, &AdjCoordSpherical);
215 MAG_Geomag(Ellip, AdjCoordSpherical, AdjCoordGeodetic, TimedMagneticModel,
216 &AdjGeoMagneticElements[0]);
217 MAG_SphericalToCartesian(AdjCoordSpherical, &x[0], &y[0], &z[0]);
218 AdjCoordGeodetic.phi = CoordGeodetic.phi - phiDelta;
219 MAG_GeodeticToSpherical(Ellip, AdjCoordGeodetic, &AdjCoordSpherical);
220 MAG_Geomag(Ellip, AdjCoordSpherical, AdjCoordGeodetic, TimedMagneticModel,
221 &AdjGeoMagneticElements[1]);
222 MAG_SphericalToCartesian(AdjCoordSpherical, &x[1], &y[1], &z[1]);
223
224 distance =
225 sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]) +
226 (z[0] - z[1]) * (z[0] - z[1]));
227 Gradient->GradPhi = MAG_GeoMagneticElementsSubtract(
228 AdjGeoMagneticElements[0], AdjGeoMagneticElements[1]);
229 Gradient->GradPhi =
230 MAG_GeoMagneticElementsScale(Gradient->GradPhi, 1 / distance);
231 AdjCoordGeodetic = MAG_CoordGeodeticAssign(CoordGeodetic);
232
233 /*Gradient along y*/
234
235 /*It is perhaps noticeable that the method here for calculation is
236 substantially different than that for the gradient along x. As we near the
237 North pole the longitude lines approach each other, and the calculation that
238 works well for latitude lines becomes unstable when 0.01 degrees represents
239 sufficiently small numbers, and fails to function correctly at all at the
240 North Pole*/
241
242 MAG_GeodeticToSpherical(Ellip, CoordGeodetic, &AdjCoordSpherical);
243 MAG_GradY(Ellip, AdjCoordSpherical, CoordGeodetic, TimedMagneticModel,
244 GeomagneticElements, &(Gradient->GradLambda));
245
246 /*Gradient along z*/
247 AdjCoordGeodetic.HeightAboveEllipsoid =
248 CoordGeodetic.HeightAboveEllipsoid + hDelta;
249 AdjCoordGeodetic.HeightAboveGeoid = CoordGeodetic.HeightAboveGeoid + hDelta;
250 MAG_GeodeticToSpherical(Ellip, AdjCoordGeodetic, &AdjCoordSpherical);
251 MAG_Geomag(Ellip, AdjCoordSpherical, AdjCoordGeodetic, TimedMagneticModel,
252 &AdjGeoMagneticElements[0]);
253 MAG_SphericalToCartesian(AdjCoordSpherical, &x[0], &y[0], &z[0]);
254 AdjCoordGeodetic.HeightAboveEllipsoid =
255 CoordGeodetic.HeightAboveEllipsoid - hDelta;
256 AdjCoordGeodetic.HeightAboveGeoid = CoordGeodetic.HeightAboveGeoid - hDelta;
257 MAG_GeodeticToSpherical(Ellip, AdjCoordGeodetic, &AdjCoordSpherical);
258 MAG_Geomag(Ellip, AdjCoordSpherical, AdjCoordGeodetic, TimedMagneticModel,
259 &AdjGeoMagneticElements[1]);
260 MAG_SphericalToCartesian(AdjCoordSpherical, &x[1], &y[1], &z[1]);
261
262 distance =
263 sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]) +
264 (z[0] - z[1]) * (z[0] - z[1]));
265 Gradient->GradZ = MAG_GeoMagneticElementsSubtract(AdjGeoMagneticElements[0],
266 AdjGeoMagneticElements[1]);
267 Gradient->GradZ = MAG_GeoMagneticElementsScale(Gradient->GradZ, 1 / distance);
268 AdjCoordGeodetic = MAG_CoordGeodeticAssign(CoordGeodetic);
269}
270
271int MAG_SetDefaults(MAGtype_Ellipsoid *Ellip, MAGtype_Geoid *Geoid)
272
273/*
274 Sets default values for WMM subroutines.
275
276 UPDATES : Ellip
277 Geoid
278
279 CALLS : none
280 */
281{
282 /* Sets WGS-84 parameters */
283 Ellip->a = 6378.137; /*semi-major axis of the ellipsoid in */
284 Ellip->b = 6356.7523142; /*semi-minor axis of the ellipsoid in */
285 Ellip->fla = 1 / 298.257223563; /* flattening */
286 Ellip->eps = sqrt(1 - (Ellip->b * Ellip->b) /
287 (Ellip->a * Ellip->a)); /*first eccentricity */
288 Ellip->epssq = (Ellip->eps * Ellip->eps); /*first eccentricity squared */
289 Ellip->re = 6371.2; /* Earth's radius */
290
291 /* Sets EGM-96 model file parameters */
292 Geoid->NumbGeoidCols =
293 1441; /* 360 degrees of longitude at 15 minute spacing */
294 Geoid->NumbGeoidRows =
295 721; /* 180 degrees of latitude at 15 minute spacing */
296 Geoid->NumbHeaderItems =
297 6; /* min, max lat, min, max long, lat, long spacing*/
298 Geoid->ScaleFactor = 4; /* 4 grid cells per degree at 15 minute spacing */
299 Geoid->NumbGeoidElevs = Geoid->NumbGeoidCols * Geoid->NumbGeoidRows;
300 Geoid->Geoid_Initialized =
301 0; /* Geoid will be initialized only if this is set to zero */
302 Geoid->UseGeoid = MAG_USE_GEOID;
303
304 return TRUE;
305} /*MAG_SetDefaults */
306
307int MAG_robustReadMagneticModel_Large(char *filename, char *filenameSV,
308 MAGtype_MagneticModel **MagneticModel) {
309 char line[MAXLINELENGTH],
310 ModelName[] = "Enhanced Magnetic Model"; /*Model Name must be no longer
311 than 31 characters*/
312 int n, nMax = 0, nMaxSV = 0, num_terms, a, epochlength = 5, i;
313 FILE *MODELFILE;
314 MODELFILE = fopen(filename, "r");
315 if (MODELFILE == 0) {
316 return 0;
317 }
318 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) {
319 return 0;
320 }
321 do {
322 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) break;
323 a = sscanf(line, "%d", &n);
324 if (n > nMax && (n < 99999 && a == 1 && n > 0)) nMax = n;
325 } while (n < 99999 && a == 1);
326 fclose(MODELFILE);
327 MODELFILE = fopen(filenameSV, "r");
328 if (MODELFILE == 0) {
329 return 0;
330 }
331 n = 0;
332 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) return 0;
333 do {
334 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) break;
335 a = sscanf(line, "%d", &n);
336 if (n > nMaxSV && (n < 99999 && a == 1 && n > 0)) nMaxSV = n;
337 } while (n < 99999 && a == 1);
338 fclose(MODELFILE);
339 num_terms = CALCULATE_NUMTERMS(nMax);
340 *MagneticModel = MAG_AllocateModelMemory(num_terms);
341 (*MagneticModel)->nMax = nMax;
342 (*MagneticModel)->nMaxSecVar = nMaxSV;
343 if (nMaxSV > 0) (*MagneticModel)->SecularVariationUsed = TRUE;
344 for (i = 0; i < num_terms; i++) {
345 (*MagneticModel)->Main_Field_Coeff_G[i] = 0;
346 (*MagneticModel)->Main_Field_Coeff_H[i] = 0;
347 (*MagneticModel)->Secular_Var_Coeff_G[i] = 0;
348 (*MagneticModel)->Secular_Var_Coeff_H[i] = 0;
349 }
350 MAG_readMagneticModel_Large(filename, filenameSV, *MagneticModel);
351 (*MagneticModel)->CoefficientFileEndDate =
352 (*MagneticModel)->epoch + epochlength;
353 strcpy((*MagneticModel)->ModelName, ModelName);
354 (*MagneticModel)->EditionDate = (*MagneticModel)->epoch;
355 return 1;
356} /*MAG_robustReadMagneticModel_Large*/
357
358int MAG_robustReadMagModels(char *filename,
359 MAGtype_MagneticModel *(*magneticmodels)[1]) {
360 int array_size = 1;
361 char line[MAXLINELENGTH];
362 int n, nMax = 0, num_terms, a;
363 FILE *MODELFILE;
364 MODELFILE = fopen(filename, "r");
365 if (MODELFILE == 0) {
366 return 0;
367 }
368 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) {
369 return 0;
370 }
371 if (line[0] == '%') {
372 MAG_readMagneticModel_SHDF(filename, magneticmodels);
373 } else if (array_size == 1) {
374 do {
375 if (NULL == fgets(line, MAXLINELENGTH, MODELFILE)) break;
376 a = sscanf(line, "%d", &n);
377 if (n > nMax && (n < 99999 && a == 1 && n > 0)) nMax = n;
378 } while (n < 99999 && a == 1);
379 num_terms = CALCULATE_NUMTERMS(nMax);
380 (*magneticmodels)[0] = MAG_AllocateModelMemory(num_terms);
381 (*magneticmodels)[0]->nMax = nMax;
382 (*magneticmodels)[0]->nMaxSecVar = nMax;
383 MAG_readMagneticModel(filename, (*magneticmodels)[0]);
384 (*magneticmodels)[0]->CoefficientFileEndDate =
385 (*magneticmodels)[0]->epoch + 5;
386
387 } else
388 return 0;
389 fclose(MODELFILE);
390 return 1;
391} /*MAG_robustReadMagModels*/
392
393/*End of Wrapper Functions*/
394
395/******************************************************************************
396 ********************************User Interface********************************
397 * This grouping consists of functions which interact with the directly with
398 * the user and are generally specific to the XXX_point.c, XXX_grid.c, and
399 * XXX_file.c programs. They deal with input from and output to the user.
400 ******************************************************************************/
401
402void MAG_Error(int control)
403
404/*This prints WMM errors.
405INPUT control Error look up number
406OUTPUT none
407CALLS : none
408
409 */
410{
411 switch (control) {
412 case 1:
413 printf("\nError allocating in MAG_LegendreFunctionMemory.\n");
414 break;
415 case 2:
416 printf("\nError allocating in MAG_AllocateModelMemory.\n");
417 break;
418 case 3:
419 printf("\nError allocating in MAG_InitializeGeoid\n");
420 break;
421 case 4:
422 printf("\nError in setting default values.\n");
423 break;
424 case 5:
425 printf("\nError initializing Geoid.\n");
426 break;
427 case 6:
428 printf("\nError opening WMM.COF\n.");
429 break;
430 case 7:
431 printf("\nError opening WMMSV.COF\n.");
432 break;
433 case 8:
434 printf("\nError reading Magnetic Model.\n");
435 break;
436 case 9:
437 printf("\nError printing Command Prompt introduction.\n");
438 break;
439 case 10:
440 printf(
441 "\nError converting from geodetic co-ordinates to spherical "
442 "co-ordinates.\n");
443 break;
444 case 11:
445 printf("\nError in time modifying the Magnetic model\n");
446 break;
447 case 12:
448 printf("\nError in Geomagnetic\n");
449 break;
450 case 13:
451 printf("\nError printing user data\n");
452 break;
453 case 14:
454 printf("\nError allocating in MAG_SummationSpecial\n");
455 break;
456 case 15:
457 printf("\nError allocating in MAG_SecVarSummationSpecial\n");
458 break;
459 case 16:
460 printf("\nError in opening EGM9615.BIN file\n");
461 break;
462 case 17:
463 printf(
464 "\nError: Latitude OR Longitude out of range in "
465 "MAG_GetGeoidHeight\n");
466 break;
467 case 18:
468 printf("\nError allocating in MAG_PcupHigh\n");
469 break;
470 case 19:
471 printf("\nError allocating in MAG_PcupLow\n");
472 break;
473 case 20:
474 printf("\nError opening coefficient file\n");
475 break;
476 case 21:
477 printf("\nError: UnitDepth too large\n");
478 break;
479 case 22:
480 printf("\nYour system needs Big endian version of EGM9615.BIN. \n");
481 printf(
482 "Please download this file from "
483 "http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml. \n");
484 printf("Replace the existing EGM9615.BIN file with the downloaded one\n");
485 break;
486 }
487} /*MAG_Error*/
488
489int MAG_GetUserGrid(MAGtype_CoordGeodetic *minimum,
490 MAGtype_CoordGeodetic *maximum, double *step_size,
491 double *a_step_size, double *step_time,
492 MAGtype_Date *StartDate, MAGtype_Date *EndDate,
493 int *ElementOption, int *PrintOption, char *OutputFile,
494 MAGtype_Geoid *Geoid)
495
496/* Prompts user to enter parameters to compute a grid - for use with the
497MAG_grid function Note: The user entries are not validated before here. The
498function populates the input variables & data structures.
499
500UPDATE : minimum Pointer to data structure with the following elements
501 double lambda; (longitude)
502 double phi; ( geodetic latitude)
503 double HeightAboveEllipsoid; (height above the ellipsoid (HaE) )
504 double HeightAboveGeoid;(height above the Geoid )
505
506 maximum -same as the above -MAG_USE_GEOID
507 step_size : double pointer : spatial step size, in decimal
508degrees a_step_size : double pointer : double altitude step size (km) step_time
509: double pointer : time step size (decimal years) StartDate : pointer to data
510structure with the following elements updates double DecimalYear; ( decimal
511years ) EndDate : Same as the above CALLS : none
512
513
514 */
515{
516 FILE *fileout;
517 char filename[] = "GridProgramDirective.txt";
518 char buffer[20];
519 int dummy;
520
521 printf("Please Enter Minimum Latitude (in decimal degrees):\n");
522 if (NULL == fgets(buffer, 20, stdin)) {
523 minimum->phi = 0;
524 printf("Unrecognized input default %lf used\n", minimum->phi);
525 } else {
526 sscanf(buffer, "%lf", &minimum->phi);
527 }
528 strcpy(buffer, "");
529 printf("Please Enter Maximum Latitude (in decimal degrees):\n");
530 if (NULL == fgets(buffer, 20, stdin)) {
531 maximum->phi = 0;
532 printf("Unrecognized input default %lf used\n", maximum->phi);
533 } else {
534 sscanf(buffer, "%lf", &maximum->phi);
535 }
536 strcpy(buffer, "");
537 printf("Please Enter Minimum Longitude (in decimal degrees):\n");
538 if (NULL == fgets(buffer, 20, stdin)) {
539 minimum->lambda = 0;
540 printf("Unrecognized input default %lf used\n", minimum->lambda);
541 } else {
542 sscanf(buffer, "%lf", &minimum->lambda);
543 }
544 strcpy(buffer, "");
545 printf("Please Enter Maximum Longitude (in decimal degrees):\n");
546 if (NULL == fgets(buffer, 20, stdin)) {
547 maximum->lambda = 0;
548 printf("Unrecognized input default %lf used\n", maximum->lambda);
549 } else {
550 sscanf(buffer, "%lf", &maximum->lambda);
551 }
552 strcpy(buffer, "");
553 printf("Please Enter Step Size (in decimal degrees):\n");
554 if (NULL == fgets(buffer, 20, stdin)) {
555 *step_size =
556 fmax(maximum->phi - minimum->phi, maximum->lambda - minimum->lambda);
557 printf("Unrecognized input default %lf used\n", *step_size);
558 } else {
559 sscanf(buffer, "%lf", step_size);
560 }
561 strcpy(buffer, "");
562 printf(
563 "Select height (default : above MSL) \n1. Above Mean Sea Level\n2. Above "
564 "WGS-84 Ellipsoid \n");
565 if (NULL == fgets(buffer, 20, stdin)) {
566 Geoid->UseGeoid = 1;
567 printf("Unrecognized option, height above MSL used.");
568 } else {
569 sscanf(buffer, "%d", &dummy);
570 if (dummy == 2)
571 Geoid->UseGeoid = 0;
572 else
573 Geoid->UseGeoid = 1;
574 }
575 strcpy(buffer, "");
576 if (Geoid->UseGeoid == 1) {
577 printf("Please Enter Minimum Height above MSL (in km):\n");
578 if (NULL == fgets(buffer, 20, stdin)) {
579 minimum->HeightAboveGeoid = 0;
580 printf("Unrecognized input default %lf used\n",
581 minimum->HeightAboveGeoid);
582 } else {
583 sscanf(buffer, "%lf", &minimum->HeightAboveGeoid);
584 }
585 strcpy(buffer, "");
586 printf("Please Enter Maximum Height above MSL (in km):\n");
587 if (NULL == fgets(buffer, 20, stdin)) {
588 maximum->HeightAboveGeoid = 0;
589 printf("Unrecognized input default %lf used\n",
590 maximum->HeightAboveGeoid);
591 } else {
592 sscanf(buffer, "%lf", &maximum->HeightAboveGeoid);
593 }
594 strcpy(buffer, "");
595
596 } else {
597 printf("Please Enter Minimum Height above the WGS-84 Ellipsoid (in km):\n");
598 if (NULL == fgets(buffer, 20, stdin)) {
599 minimum->HeightAboveGeoid = 0;
600 // printf("Unrecognized input default %lf used\n",
601 // minimum->HeightAboveGeoid);
602 } else {
603 sscanf(buffer, "%lf", &minimum->HeightAboveGeoid);
604 }
605 minimum->HeightAboveEllipsoid = minimum->HeightAboveGeoid;
606 strcpy(buffer, "");
607 printf("Please Enter Maximum Height above the WGS-84 Ellipsoid (in km):\n");
608 if (NULL == fgets(buffer, 20, stdin)) {
609 maximum->HeightAboveGeoid = 0;
610 // printf("Unrecognized input default %lf used\n",
611 // maximum->HeightAboveGeoid);
612 } else {
613 sscanf(buffer, "%lf", &maximum->HeightAboveGeoid);
614 }
615 maximum->HeightAboveEllipsoid = maximum->HeightAboveGeoid;
616 strcpy(buffer, "");
617 }
618 printf("Please Enter height step size (in km):\n");
619 if (NULL == fgets(buffer, 20, stdin)) {
620 *a_step_size = maximum->HeightAboveGeoid - minimum->HeightAboveGeoid;
621 printf("Unrecognized input default %lf used\n", *a_step_size);
622 } else {
623 sscanf(buffer, "%lf", a_step_size);
624 }
625 strcpy(buffer, "");
626 printf("\nPlease Enter the decimal year starting time:\n");
627 while (NULL == fgets(buffer, 20, stdin)) {
628 printf("\nUnrecognized input, please re-enter a decimal year\n");
629 }
630 sscanf(buffer, "%lf", &StartDate->DecimalYear);
631 strcpy(buffer, "");
632 printf("Please Enter the decimal year ending time:\n");
633 while (NULL == fgets(buffer, 20, stdin)) {
634 printf("\nUnrecognized input, please re-enter a decimal year\n");
635 }
636 sscanf(buffer, "%lf", &EndDate->DecimalYear);
637 strcpy(buffer, "");
638 printf("Please Enter the time step size:\n");
639 if (NULL == fgets(buffer, 20, stdin)) {
640 *step_time = EndDate->DecimalYear - StartDate->DecimalYear;
641 printf("Unrecognized input, default of %lf used\n", *step_time);
642 } else {
643 sscanf(buffer, "%lf", step_time);
644 }
645 strcpy(buffer, "");
646 printf("Enter a geomagnetic element to print. Your options are:\n");
647 printf(
648 " 1. Declination 9. Ddot\n 2. Inclination 10. Idot\n 3. F 11. "
649 "Fdot\n 4. H 12. Hdot\n 5. X 13. Xdot\n 6. Y 14. Ydot\n 7. Z "
650 " 15. Zdot\n 8. GV 16. GVdot\nFor gradients enter: 17\n");
651 if (NULL == fgets(buffer, 20, stdin)) {
652 *ElementOption = 1;
653 printf("Unrecognized input, default of %d used\n", *ElementOption);
654 }
655 sscanf(buffer, "%d", ElementOption);
656 strcpy(buffer, "");
657 if (*ElementOption == 17) {
658 printf("Enter a gradient element to print. Your options are:\n");
659 printf(" 1. dX/dphi \t2. dY/dphi \t3. dZ/dphi\n");
660 printf(" 4. dX/dlambda \t5. dY/dlambda \t6. dZ/dlambda\n");
661 printf(" 7. dX/dz \t8. dY/dz \t9. dZ/dz\n");
662 strcpy(buffer, "");
663 if (NULL == fgets(buffer, 20, stdin)) {
664 *ElementOption = 1;
665 printf("Unrecognized input, default of %d used\n", *ElementOption);
666 } else {
667 sscanf(buffer, "%d", ElementOption);
668 }
669 strcpy(buffer, "");
670 *ElementOption += 16;
671 }
672 printf("Select output :\n");
673 printf(" 1. Print to a file \n 2. Print to Screen\n");
674 if (NULL == fgets(buffer, 20, stdin)) {
675 *PrintOption = 2;
676 printf("Unrecognized input, default of printing to screen\n");
677 } else {
678 sscanf(buffer, "%d", PrintOption);
679 }
680 strcpy(buffer, "");
681 fileout = fopen(filename, "a");
682 if (*PrintOption == 1) {
683 printf(
684 "Please enter output filename\nfor default ('GridResults.txt') press "
685 "enter:\n");
686 if (NULL == fgets(buffer, 20, stdin) || strlen(buffer) <= 1) {
687 strcpy(OutputFile, "GridResults.txt");
688 fprintf(fileout, "\nResults printed in: GridResults.txt\n");
689 strcpy(OutputFile, "GridResults.txt");
690 } else {
691 sscanf(buffer, "%s", OutputFile);
692 fprintf(fileout, "\nResults printed in: %s\n", OutputFile);
693 }
694 /*strcpy(OutputFile, buffer);*/
695 strcpy(buffer, "");
696 /*sscanf(buffer, "%s", OutputFile);*/
697 } else
698 fprintf(fileout, "\nResults printed in Console\n");
699 fprintf(
700 fileout,
701 "Minimum Latitude: %f\t\tMaximum Latitude: %f\t\tStep Size: %f\nMinimum "
702 "Longitude: %f\t\tMaximum Longitude: %f\t\tStep Size: %f\n",
703 minimum->phi, maximum->phi, *step_size, minimum->lambda, maximum->lambda,
704 *step_size);
705 if (Geoid->UseGeoid == 1)
706 fprintf(fileout,
707 "Minimum Altitude above MSL: %f\tMaximum Altitude above MSL: "
708 "%f\tStep Size: %f\n",
709 minimum->HeightAboveGeoid, maximum->HeightAboveGeoid, *a_step_size);
710 else
711 fprintf(fileout,
712 "Minimum Altitude above WGS-84 Ellipsoid: %f\tMaximum Altitude "
713 "above WGS-84 Ellipsoid: %f\tStep Size: %f\n",
714 minimum->HeightAboveEllipsoid, maximum->HeightAboveEllipsoid,
715 *a_step_size);
716 fprintf(fileout,
717 "Starting Date: %f\t\tEnding Date: %f\t\tStep Time: %f\n\n\n",
718 StartDate->DecimalYear, EndDate->DecimalYear, *step_time);
719 fclose(fileout);
720 return TRUE;
721}
722
723#ifndef OPENCPN
724int MAG_GetUserInput(MAGtype_MagneticModel *MagneticModel, MAGtype_Geoid *Geoid,
725 MAGtype_CoordGeodetic *CoordGeodetic,
726 MAGtype_Date *MagneticDate)
727
728/*
729This prompts the user for coordinates, and accepts many entry formats.
730It takes the MagneticModel and Geoid as input and outputs the Geographic
731coordinates and Date as objects. Returns 0 when the user wants to exit and 1 if
732the user enters valid input data. INPUT : MagneticModel : Data structure with
733the following elements used here double epoch; Base time of Geomagnetic
734model epoch (yrs) : Geoid Pointer to data structure MAGtype_Geoid (used for
735converting HeightAboveGeoid to HeightABoveEllipsoid
736
737OUTPUT: CoordGeodetic : Pointer to data structure. Following elements are
738updated double lambda; (longitude) double phi; ( geodetic latitude) double
739HeightAboveEllipsoid; (height above the ellipsoid (HaE) ) double
740HeightAboveGeoid;(height above the Geoid )
741
742 MagneticDate : Pointer to data structure MAGtype_Date with the
743following elements updated int Year; (If user directly enters decimal year this
744field is not populated) int Month;(If user directly enters decimal year this
745field is not populated) int Day; (If user directly enters decimal year this
746field is not populated) double DecimalYear; decimal years
747
748CALLS: MAG_DMSstringToDegree(buffer, &CoordGeodetic->lambda); (The program uses
749this to convert the string into a decimal longitude.)
750 MAG_ValidateDMSstringlong(buffer, Error_Message)
751 MAG_ValidateDMSstringlat(buffer, Error_Message)
752 MAG_Warnings
753 MAG_ConvertGeoidToEllipsoidHeight
754 MAG_DateToYear
755
756 */
757{
758 char Error_Message[255];
759 char buffer[40];
760 int i, j, a, b, c, done = 0;
761 double lat_bound[2] = {LAT_BOUND_MIN, LAT_BOUND_MAX};
762 double lon_bound[2] = {LON_BOUND_MIN, LON_BOUND_MAX};
763 int alt_bound[2] = {ALT_BOUND_MIN, NO_ALT_MAX};
764 char *Qstring = malloc(sizeof(char) * 1028);
765 strcpy(buffer, ""); /*Clear the input */
766 strcpy(Qstring,
767 "\nPlease enter latitude\nNorth latitude positive, For example:\n30, "
768 "30, 30 (D,M,S) or 30.508 (Decimal Degrees) (both are north)\n");
769 MAG_GetDeg(Qstring, &CoordGeodetic->phi, lat_bound);
770 strcpy(buffer, ""); /*Clear the input*/
771 strcpy(Qstring,
772 "\nPlease enter longitude\nEast longitude positive, West negative. "
773 "For example:\n-100.5 or -100, 30, 0 for 100.5 degrees west\n");
774 MAG_GetDeg(Qstring, &CoordGeodetic->lambda, lon_bound);
775
776 strcpy(Qstring,
777 "\nPlease enter height above mean sea level (in kilometers):\n[For "
778 "height above WGS-84 ellipsoid prefix E, for example (E20.1)]\n");
779 if (MAG_GetAltitude(Qstring, Geoid, CoordGeodetic, alt_bound, FALSE) ==
780 USER_GAVE_UP)
781 return FALSE;
782 strcpy(buffer, "");
783 printf(
784 "\nPlease enter the decimal year or calendar date\n (YYYY.yyy, MM DD "
785 "YYYY or MM/DD/YYYY):\n");
786 while (NULL == fgets(buffer, 40, stdin)) {
787 printf(
788 "\nPlease enter the decimal year or calendar date\n (YYYY.yyy, MM DD "
789 "YYYY or MM/DD/YYYY):\n");
790 }
791 for (i = 0, done = 0; i <= 40 && !done; i++) {
792 if (buffer[i] == '.') {
793 j = sscanf(buffer, "%lf", &MagneticDate->DecimalYear);
794 if (j == 1)
795 done = 1;
796 else
797 buffer[i] = '\0';
798 }
799 if (buffer[i] == '/') {
800 sscanf(buffer, "%d/%d/%d", &MagneticDate->Month, &MagneticDate->Day,
801 &MagneticDate->Year);
802 if (!MAG_DateToYear(MagneticDate, Error_Message)) {
803 printf("%s", Error_Message);
804 printf(
805 "\nPlease re-enter Date in MM/DD/YYYY or MM DD YYYY format, or as "
806 "a decimal year\n");
807 while (NULL == fgets(buffer, 40, stdin)) {
808 printf(
809 "\nPlease re-enter Date in MM/DD/YYYY or MM DD YYYY format, or "
810 "as a decimal year\n");
811 }
812 i = 0;
813 } else
814 done = 1;
815 }
816 if ((buffer[i] == ' ' && buffer[i + 1] != '/') || buffer[i] == '\0') {
817 if (3 == sscanf(buffer, "%d %d %d", &a, &b, &c)) {
818 MagneticDate->Month = a;
819 MagneticDate->Day = b;
820 MagneticDate->Year = c;
821 MagneticDate->DecimalYear = 99999;
822 } else if (1 == sscanf(buffer, "%d %d %d", &a, &b, &c)) {
823 MagneticDate->DecimalYear = a;
824 done = 1;
825 }
826 if (!(MagneticDate->DecimalYear == a)) {
827 if (!MAG_DateToYear(MagneticDate, Error_Message)) {
828 printf("%s", Error_Message);
829 strcpy(buffer, "");
830 printf(
831 "\nError encountered, please re-enter Date in MM/DD/YYYY or MM "
832 "DD YYYY format, or as a decimal year\n");
833 while (NULL == fgets(buffer, 40, stdin)) {
834 printf(
835 "\nError encountered, please re-enter Date in MM/DD/YYYY or MM "
836 "DD YYYY format, or as a decimal year\n");
837 }
838 i = -1;
839 } else
840 done = 1;
841 }
842 }
843 if (buffer[i] == '\0' && i != -1 && done != 1) {
844 strcpy(buffer, "");
845 printf(
846 "\nError encountered, please re-enter as MM/DD/YYYY, MM DD YYYY, or "
847 "as YYYY.yyy:\n");
848 while (NULL == fgets(buffer, 40, stdin)) {
849 printf(
850 "\nError encountered, please re-enter as MM/DD/YYYY, MM DD YYYY, "
851 "or as YYYY.yyy:\n");
852 }
853 i = -1;
854 }
855 if (done) {
856 if (MagneticDate->DecimalYear > MagneticModel->CoefficientFileEndDate ||
857 MagneticDate->DecimalYear < MagneticModel->epoch) {
858 switch (MAG_Warnings(4, MagneticDate->DecimalYear, MagneticModel)) {
859 case 0:
860 return 0;
861 case 1:
862 done = 0;
863 i = -1;
864 strcpy(buffer, "");
865 printf(
866 "\nPlease enter the decimal year or calendar date\n (YYYY.yyy, "
867 "MM DD YYYY or MM/DD/YYYY):\n");
868 while (NULL == fgets(buffer, 40, stdin)) {
869 printf(
870 "\nPlease enter the decimal year or calendar date\n "
871 "(YYYY.yyy, MM DD YYYY or MM/DD/YYYY):\n");
872 }
873 break;
874 case 2:
875 break;
876 }
877 }
878 }
879 }
880 free(Qstring);
881 return TRUE;
882} /*MAG_GetUserInput*/
883
884#endif /* OPENCPN */
885
886void MAG_PrintGradient(MAGtype_Gradient Gradient) {
887 printf("\nGradient\n");
888 printf("\n Northward Eastward Downward\n");
889 printf("X: %7.1f nT/km %9.1f nT/km %9.1f nT/km \n",
890 Gradient.GradPhi.X, Gradient.GradLambda.X, Gradient.GradZ.X);
891 printf("Y: %7.1f nT/km %9.1f nT/km %9.1f nT/km \n",
892 Gradient.GradPhi.Y, Gradient.GradLambda.Y, Gradient.GradZ.Y);
893 printf("Z: %7.1f nT/km %9.1f nT/km %9.1f nT/km \n",
894 Gradient.GradPhi.Z, Gradient.GradLambda.Z, Gradient.GradZ.Z);
895 printf("H: %7.1f nT/km %9.1f nT/km %9.1f nT/km \n",
896 Gradient.GradPhi.H, Gradient.GradLambda.H, Gradient.GradZ.H);
897 printf("F: %7.1f nT/km %9.1f nT/km %9.1f nT/km \n",
898 Gradient.GradPhi.F, Gradient.GradLambda.F, Gradient.GradZ.F);
899 printf("Declination: %7.2f min/km %8.2f min/km %8.2f min/km \n",
900 Gradient.GradPhi.Decl * 60, Gradient.GradLambda.Decl * 60,
901 Gradient.GradZ.Decl * 60);
902 printf("Inclination: %7.2f min/km %8.2f min/km %8.2f min/km \n",
903 Gradient.GradPhi.Incl * 60, Gradient.GradLambda.Incl * 60,
904 Gradient.GradZ.Incl * 60);
905}
906
907void MAG_PrintUserData(MAGtype_GeoMagneticElements GeomagElements,
908 MAGtype_CoordGeodetic SpaceInput, MAGtype_Date TimeInput,
909 MAGtype_MagneticModel *MagneticModel,
910 MAGtype_Geoid *Geoid)
911/* This function prints the results in Geomagnetic Elements for a point
912calculation. It takes the calculated
913 * Geomagnetic elements "GeomagElements" as input.
914 * As well as the coordinates, date, and Magnetic Model.
915INPUT : GeomagElements : Data structure MAGtype_GeoMagneticElements with the
916following elements double Decl; (Angle between the magnetic field vector and
917true north, positive east) double Incl; Angle between the magnetic field vector
918and the horizontal plane, positive down double F; Magnetic Field Strength double
919H; Horizontal Magnetic Field Strength double X; Northern component of the
920magnetic field vector double Y; Eastern component of the magnetic field vector
921 double Z; Downward component of the magnetic field
922vector4 double Decldot; Yearly Rate of change in declination double Incldot;
923Yearly Rate of change in inclination double Fdot; Yearly rate of change in
924Magnetic field strength double Hdot; Yearly rate of change in horizontal field
925strength double Xdot; Yearly rate of change in the northern component double
926Ydot; Yearly rate of change in the eastern component double Zdot; Yearly rate of
927change in the downward component double GVdot;Yearly rate of chnage in grid
928variation CoordGeodetic Pointer to the data structure with the following
929elements double lambda; (longitude) double phi; ( geodetic latitude) double
930HeightAboveEllipsoid; (height above the ellipsoid (HaE) ) double
931HeightAboveGeoid;(height above the Geoid ) TimeInput : data structure
932MAGtype_Date with the following elements int Year; int Month; int Day; double
933DecimalYear; decimal years MagneticModel : data structure with the
934following elements double EditionDate; double epoch; Base time of
935Geomagnetic model epoch (yrs) char ModelName[20]; double *Main_Field_Coeff_G;
936C - Gauss coefficients of main geomagnetic model (nT) double
937*Main_Field_Coeff_H; C - Gauss coefficients of main geomagnetic model
938(nT) double *Secular_Var_Coeff_G; CD - Gauss coefficients of secular
939geomagnetic model (nT/yr) double *Secular_Var_Coeff_H; CD - Gauss coefficients
940of secular geomagnetic model (nT/yr) int nMax; Maximum degree of spherical
941harmonic model int nMaxSecVar; Maxumum degree of spherical harmonic secular
942model int SecularVariationUsed; Whether or not the magnetic secular variation
943vector will be needed by program OUTPUT : none
944 */
945{
946 char DeclString[100];
947 char InclString[100];
948 MAG_DegreeToDMSstring(GeomagElements.Incl, 2, InclString);
949 if (GeomagElements.H < 6000 && GeomagElements.H > 2000)
950 MAG_Warnings(1, GeomagElements.H, MagneticModel);
951 if (GeomagElements.H < 2000) MAG_Warnings(2, GeomagElements.H, MagneticModel);
952 if (MagneticModel->SecularVariationUsed == TRUE) {
953 MAG_DegreeToDMSstring(GeomagElements.Decl, 2, DeclString);
954 printf("\n Results For \n\n");
955 if (SpaceInput.phi < 0)
956 printf("Latitude %.2fS\n", -SpaceInput.phi);
957 else
958 printf("Latitude %.2fN\n", SpaceInput.phi);
959 if (SpaceInput.lambda < 0)
960 printf("Longitude %.2fW\n", -SpaceInput.lambda);
961 else
962 printf("Longitude %.2fE\n", SpaceInput.lambda);
963 if (Geoid->UseGeoid == 1)
964 printf("Altitude: %.2f Kilometers above mean sea level\n",
965 SpaceInput.HeightAboveGeoid);
966 else
967 printf("Altitude: %.2f Kilometers above the WGS-84 ellipsoid\n",
968 SpaceInput.HeightAboveEllipsoid);
969 printf("Date: %.1f\n", TimeInput.DecimalYear);
970 printf("\n Main Field\t\t\tSecular Change\n");
971 printf("F = %-9.1f nT\t\t Fdot = %.1f\tnT/yr\n", GeomagElements.F,
972 GeomagElements.Fdot);
973 printf("H = %-9.1f nT\t\t Hdot = %.1f\tnT/yr\n", GeomagElements.H,
974 GeomagElements.Hdot);
975 printf("X = %-9.1f nT\t\t Xdot = %.1f\tnT/yr\n", GeomagElements.X,
976 GeomagElements.Xdot);
977 printf("Y = %-9.1f nT\t\t Ydot = %.1f\tnT/yr\n", GeomagElements.Y,
978 GeomagElements.Ydot);
979 printf("Z = %-9.1f nT\t\t Zdot = %.1f\tnT/yr\n", GeomagElements.Z,
980 GeomagElements.Zdot);
981 if (GeomagElements.Decl < 0)
982 printf("Decl =%20s (WEST)\t Ddot = %.1f\tMin/yr\n", DeclString,
983 60 * GeomagElements.Decldot);
984 else
985 printf("Decl =%20s (EAST)\t Ddot = %.1f\tMin/yr\n", DeclString,
986 60 * GeomagElements.Decldot);
987 if (GeomagElements.Incl < 0)
988 printf("Incl =%20s (UP)\t Idot = %.1f\tMin/yr\n", InclString,
989 60 * GeomagElements.Incldot);
990 else
991 printf("Incl =%20s (DOWN)\t Idot = %.1f\tMin/yr\n", InclString,
992 60 * GeomagElements.Incldot);
993 } else {
994 MAG_DegreeToDMSstring(GeomagElements.Decl, 2, DeclString);
995 printf("\n Results For \n\n");
996 if (SpaceInput.phi < 0)
997 printf("Latitude %.2fS\n", -SpaceInput.phi);
998 else
999 printf("Latitude %.2fN\n", SpaceInput.phi);
1000 if (SpaceInput.lambda < 0)
1001 printf("Longitude %.2fW\n", -SpaceInput.lambda);
1002 else
1003 printf("Longitude %.2fE\n", SpaceInput.lambda);
1004 if (Geoid->UseGeoid == 1)
1005 printf("Altitude: %.2f Kilometers above MSL\n",
1006 SpaceInput.HeightAboveGeoid);
1007 else
1008 printf("Altitude: %.2f Kilometers above WGS-84 Ellipsoid\n",
1009 SpaceInput.HeightAboveEllipsoid);
1010 printf("Date: %.1f\n", TimeInput.DecimalYear);
1011 printf("\n Main Field\n");
1012 printf("F = %-9.1f nT\n", GeomagElements.F);
1013 printf("H = %-9.1f nT\n", GeomagElements.H);
1014 printf("X = %-9.1f nT\n", GeomagElements.X);
1015 printf("Y = %-9.1f nT\n", GeomagElements.Y);
1016 printf("Z = %-9.1f nT\n", GeomagElements.Z);
1017 if (GeomagElements.Decl < 0)
1018 printf("Decl =%20s (WEST)\n", DeclString);
1019 else
1020 printf("Decl =%20s (EAST)\n", DeclString);
1021 if (GeomagElements.Incl < 0)
1022 printf("Incl =%20s (UP)\n", InclString);
1023 else
1024 printf("Incl =%20s (DOWN)\n", InclString);
1025 }
1026
1027 if (SpaceInput.phi <= -55 || SpaceInput.phi >= 55)
1028 /* Print Grid Variation */
1029 {
1030 MAG_DegreeToDMSstring(GeomagElements.GV, 2, InclString);
1031 printf("\n\n Grid variation =%20s\n", InclString);
1032 }
1033
1034} /*MAG_PrintUserData*/
1035
1036int MAG_ValidateDMSstring(char *input, int min, int max, char *Error)
1037
1038/* Validates a latitude DMS string, and returns 1 for a success and returns 0
1039for a failure. It copies an error message to the Error string in the event of a
1040failure.
1041
1042INPUT : input (DMS string)
1043OUTPUT : Error : Error string
1044CALLS : none
1045 */
1046{
1047 int degree, minute, second, j = 0, n, max_minute = 60, max_second = 60;
1048 int i;
1049 degree = -1000;
1050 minute = -1;
1051 second = -1;
1052 n = (int)strlen(input);
1053
1054 for (i = 0; i <= n - 1; i++) /*tests for legal characters*/
1055 {
1056 if ((input[i] < '0' || input[i] > '9') &&
1057 (input[i] != ',' && input[i] != ' ' && input[i] != '-' &&
1058 input[i] != '\0' && input[i] != '\n')) {
1059 strcpy(Error,
1060 "\nError: Input contains an illegal character, legal characters "
1061 "for Degree, Minute, Second format are:\n '0-9' ',' '-' '[space]' "
1062 "'[Enter]'\n");
1063 return FALSE;
1064 }
1065 if (input[i] == ',') j++;
1066 }
1067 if (j == 2)
1068 j = sscanf(input, "%d, %d, %d", &degree, &minute,
1069 &second); /*tests for legal formatting and range*/
1070 else
1071 j = sscanf(input, "%d %d %d", &degree, &minute, &second);
1072 if (j == 1) {
1073 minute = 0;
1074 second = 0;
1075 j = 3;
1076 }
1077 if (j != 3) {
1078 strcpy(Error,
1079 "\nError: Not enough numbers used for Degrees, Minutes, Seconds "
1080 "format\n or they were incorrectly formatted\n The legal format is "
1081 "DD,MM,SS or DD MM SS\n");
1082 return FALSE;
1083 }
1084 if (degree > max || degree < min) {
1085 sprintf(Error,
1086 "\nError: Degree input is outside legal range\n The legal range is "
1087 "from %d to %d\n",
1088 min, max);
1089 return FALSE;
1090 }
1091 if (degree == max || degree == min) max_minute = 0;
1092 if (minute > max_minute || minute < 0) {
1093 strcpy(Error,
1094 "\nError: Minute input is outside legal range\n The legal minute "
1095 "range is from 0 to 60\n");
1096 return FALSE;
1097 }
1098 if (minute == max_minute) max_second = 0;
1099 if (second > max_second || second < 0) {
1100 strcpy(Error,
1101 "\nError: Second input is outside legal range\n The legal second "
1102 "range is from 0 to 60\n");
1103 return FALSE;
1104 }
1105 return TRUE;
1106} /*MAG_ValidateDMSstring*/
1107
1108int MAG_Warnings(int control, double value,
1109 MAGtype_MagneticModel *MagneticModel)
1110
1111/*Return value 0 means end program, Return value 1 means get new data, Return
1112value 2 means continue. This prints a warning to the screen determined by the
1113control integer. It also takes the value of the parameter causing the warning as
1114a double. This is unnecessary for some warnings. It requires the MagneticModel
1115to determine the current epoch.
1116
1117 INPUT control :int : (Warning number)
1118 value : double: Magnetic field strength
1119 MagneticModel
1120OUTPUT : none
1121CALLS : none
1122
1123 */
1124{
1125 char ans[20];
1126 strcpy(ans, "");
1127
1128 switch (control) {
1129 case 1: /* Horizontal Field strength low */
1130 do {
1131 printf(
1132 "\nCaution: location is approaching the blackout zone around the "
1133 "magnetic pole as\n");
1134 printf(" defined by the WMM military specification \n");
1135 printf(
1136 " "
1137 "(https://www.ngdc.noaa.gov/geomag/WMM/data/MIL-PRF-89500B.pdf). "
1138 "Compass\n");
1139 printf(" accuracy may be degraded in this region.\n");
1140 printf("Press enter to continue...\n");
1141 } while (NULL == fgets(ans, 20, stdin));
1142 break;
1143 case 2: /* Horizontal Field strength very low */
1144 do {
1145 printf(
1146 "\nWarning: location is in the blackout zone around the magnetic "
1147 "pole as defined\n");
1148 printf(" by the WMM military specification \n");
1149 printf(
1150 " "
1151 "(https://www.ngdc.noaa.gov/geomag/WMM/data/MIL-PRF-89500B.pdf). "
1152 "Compass\n");
1153 printf(" accuracy is highly degraded in this region.\n");
1154 } while (NULL == fgets(ans, 20, stdin));
1155 break;
1156 case 3: /* Elevation outside the recommended range */
1157 printf(
1158 "\nWarning: The value you have entered of %.1f km for the elevation "
1159 "is outside of the recommended range.\n Elevations above -10.0 km "
1160 "are recommended for accurate results. \n",
1161 value);
1162 while (1) {
1163 printf(
1164 "\nPlease press 'C' to continue, 'G' to get new data or 'X' to "
1165 "exit...\n");
1166 while (NULL == fgets(ans, 20, stdin)) {
1167 printf("\nInvalid input\n");
1168 }
1169 switch (ans[0]) {
1170 case 'X':
1171 case 'x':
1172 return 0;
1173 case 'G':
1174 case 'g':
1175 return 1;
1176 case 'C':
1177 case 'c':
1178 return 2;
1179 default:
1180 printf("\nInvalid input %c\n", ans[0]);
1181 break;
1182 }
1183 }
1184 break;
1185
1186 case 4: /*Date outside the recommended range*/
1187 printf(
1188 "\nWARNING - TIME EXTENDS BEYOND INTENDED USAGE RANGE\n CONTACT NCEI "
1189 "FOR PRODUCT UPDATES:\n");
1190 printf(" National Centers for Environmental Information\n");
1191 printf(" NOAA E/NE42\n");
1192 printf(" 325 Broadway\n");
1193 printf("\n Boulder, CO 80305 USA");
1194 printf(" Attn: Manoj Nair or Arnaud Chulliat\n");
1195 printf(" Phone: (303) 497-4642 or -6522\n");
1196 printf(" Email: geomag.models@noaa.gov\n");
1197 printf(" Web: http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml\n");
1198 printf("\n VALID RANGE = %d - %d\n", (int)MagneticModel->epoch,
1199 (int)MagneticModel->CoefficientFileEndDate);
1200 printf(" TIME = %f\n", value);
1201 while (1) {
1202 printf(
1203 "\nPlease press 'C' to continue, 'N' to enter new data or 'X' to "
1204 "exit...\n");
1205 while (NULL == fgets(ans, 20, stdin)) {
1206 printf("\nInvalid input\n");
1207 }
1208 switch (ans[0]) {
1209 case 'X':
1210 case 'x':
1211 return 0;
1212 case 'N':
1213 case 'n':
1214 return 1;
1215 case 'C':
1216 case 'c':
1217 return 2;
1218 default:
1219 printf("\nInvalid input %c\n", ans[0]);
1220 break;
1221 }
1222 }
1223 break;
1224 case 5: /*Elevation outside the allowable range*/
1225 printf(
1226 "\nError: The value you have entered of %f km for the elevation is "
1227 "outside of the recommended range.\n Elevations above -10.0 km are "
1228 "recommended for accurate results. \n",
1229 value);
1230 while (1) {
1231 printf(
1232 "\nPlease press 'C' to continue, 'G' to get new data or 'X' to "
1233 "exit...\n");
1234 while (NULL == fgets(ans, 20, stdin)) {
1235 printf("\nInvalid input\n");
1236 }
1237 switch (ans[0]) {
1238 case 'X':
1239 case 'x':
1240 return 0;
1241 case 'G':
1242 case 'g':
1243 return 1;
1244 case 'C':
1245 case 'c':
1246 return 2;
1247 default:
1248 printf("\nInvalid input %c\n", ans[0]);
1249 break;
1250 }
1251 }
1252 break;
1253 }
1254 return 2;
1255} /*MAG_Warnings*/
1256
1257/*End of User Interface functions*/
1258
1259/******************************************************************************
1260 ********************************Memory and File Processing********************
1261 * This grouping consists of functions that read coefficient files into the
1262 * memory, allocate memory, free memory or print models into coefficient files.
1263 ******************************************************************************/
1264
1265MAGtype_LegendreFunction *MAG_AllocateLegendreFunctionMemory(int NumTerms)
1266
1267/* Allocate memory for Associated Legendre Function data types.
1268 Should be called before computing Associated Legendre Functions.
1269
1270 INPUT: NumTerms : int : Total number of spherical harmonic coefficients in the
1271model
1272
1273
1274 OUTPUT: Pointer to data structure MAGtype_LegendreFunction with the
1275following elements double *Pcup; ( pointer to store Legendre Function )
1276 double *dPcup; ( pointer to store Derivative of
1277Legendre function )
1278
1279 FALSE: Failed to allocate memory
1280
1281CALLS : none
1282
1283 */
1284{
1285 MAGtype_LegendreFunction *LegendreFunction;
1286
1287 LegendreFunction =
1288 (MAGtype_LegendreFunction *)calloc(1, sizeof(MAGtype_LegendreFunction));
1289
1290 if (!LegendreFunction) {
1291 MAG_Error(1);
1292 return NULL;
1293 }
1294 LegendreFunction->Pcup = (double *)malloc((NumTerms + 1) * sizeof(double));
1295 if (LegendreFunction->Pcup == 0) {
1296 MAG_Error(1);
1297 return NULL;
1298 }
1299 LegendreFunction->dPcup = (double *)malloc((NumTerms + 1) * sizeof(double));
1300 if (LegendreFunction->dPcup == 0) {
1301 MAG_Error(1);
1302 return NULL;
1303 }
1304 return LegendreFunction;
1305} /*MAGtype_LegendreFunction*/
1306
1307MAGtype_MagneticModel *MAG_AllocateModelMemory(int NumTerms)
1308
1309/* Allocate memory for WMM Coefficients
1310 * Should be called before reading the model file *
1311
1312 INPUT: NumTerms : int : Total number of spherical harmonic coefficients in the
1313model
1314
1315
1316 OUTPUT: Pointer to data structure MAGtype_MagneticModel with the following
1317elements double EditionDate; double epoch; Base time of Geomagnetic model
1318epoch (yrs) char ModelName[20]; double *Main_Field_Coeff_G; C - Gauss
1319coefficients of main geomagnetic model (nT) double *Main_Field_Coeff_H; C -
1320Gauss coefficients of main geomagnetic model (nT) double *Secular_Var_Coeff_G;
1321CD - Gauss coefficients of secular geomagnetic model (nT/yr) double
1322*Secular_Var_Coeff_H; CD - Gauss coefficients of secular geomagnetic model
1323(nT/yr) int nMax; Maximum degree of spherical harmonic model int nMaxSecVar;
1324Maxumum degree of spherical harmonic secular model int SecularVariationUsed;
1325Whether or not the magnetic secular variation vector will be needed by program
1326
1327 FALSE: Failed to allocate memory
1328CALLS : none
1329 */
1330{
1331 MAGtype_MagneticModel *MagneticModel;
1332 int i;
1333
1334 MagneticModel =
1335 (MAGtype_MagneticModel *)calloc(1, sizeof(MAGtype_MagneticModel));
1336
1337 if (MagneticModel == NULL) {
1338 MAG_Error(2);
1339 return NULL;
1340 }
1341
1342 MagneticModel->Main_Field_Coeff_G =
1343 (double *)malloc((NumTerms + 1) * sizeof(double));
1344
1345 if (MagneticModel->Main_Field_Coeff_G == NULL) {
1346 MAG_Error(2);
1347 return NULL;
1348 }
1349
1350 MagneticModel->Main_Field_Coeff_H =
1351 (double *)malloc((NumTerms + 1) * sizeof(double));
1352
1353 if (MagneticModel->Main_Field_Coeff_H == NULL) {
1354 MAG_Error(2);
1355 return NULL;
1356 }
1357 MagneticModel->Secular_Var_Coeff_G =
1358 (double *)malloc((NumTerms + 1) * sizeof(double));
1359 if (MagneticModel->Secular_Var_Coeff_G == NULL) {
1360 MAG_Error(2);
1361 return NULL;
1362 }
1363 MagneticModel->Secular_Var_Coeff_H =
1364 (double *)malloc((NumTerms + 1) * sizeof(double));
1365 if (MagneticModel->Secular_Var_Coeff_H == NULL) {
1366 MAG_Error(2);
1367 return NULL;
1368 }
1369 MagneticModel->CoefficientFileEndDate = 0;
1370 MagneticModel->EditionDate = 0;
1371 strcpy(MagneticModel->ModelName, "");
1372 MagneticModel->SecularVariationUsed = 0;
1373 MagneticModel->epoch = 0;
1374 MagneticModel->nMax = 0;
1375 MagneticModel->nMaxSecVar = 0;
1376
1377 for (i = 0; i < NumTerms; i++) {
1378 MagneticModel->Main_Field_Coeff_G[i] = 0;
1379 MagneticModel->Main_Field_Coeff_H[i] = 0;
1380 MagneticModel->Secular_Var_Coeff_G[i] = 0;
1381 MagneticModel->Secular_Var_Coeff_H[i] = 0;
1382 }
1383
1384 return MagneticModel;
1385
1386} /*MAG_AllocateModelMemory*/
1387
1388MAGtype_SphericalHarmonicVariables *MAG_AllocateSphVarMemory(int nMax) {
1389 MAGtype_SphericalHarmonicVariables *SphVariables;
1390 SphVariables = (MAGtype_SphericalHarmonicVariables *)calloc(
1391 1, sizeof(MAGtype_SphericalHarmonicVariables));
1392 SphVariables->RelativeRadiusPower =
1393 (double *)malloc((nMax + 1) * sizeof(double));
1394 SphVariables->cos_mlambda = (double *)malloc((nMax + 1) * sizeof(double));
1395 SphVariables->sin_mlambda = (double *)malloc((nMax + 1) * sizeof(double));
1396 return SphVariables;
1397} /*MAG_AllocateSphVarMemory*/
1398
1399void MAG_AssignHeaderValues(MAGtype_MagneticModel *model,
1400 char values[][MAXLINELENGTH]) {
1401 /* MAGtype_Date releasedate; */
1402 strcpy(model->ModelName, values[MODELNAME]);
1403 /* releasedate.Year = 0;
1404 releasedate.Day = 0;
1405 releasedate.Month = 0;
1406 releasedate.DecimalYear = 0;
1407 sscanf(values[RELEASEDATE],"%d-%d-%d",&releasedate.Year,&releasedate.Month,&releasedate.Day);
1408 if(MAG_DateToYear (&releasedate, NULL))
1409 model->EditionDate = releasedate.DecimalYear;*/
1410 model->epoch = atof(values[MODELSTARTYEAR]);
1411 model->nMax = atoi(values[INTSTATICDEG]);
1412 model->nMaxSecVar = atoi(values[INTSECVARDEG]);
1413 model->CoefficientFileEndDate = atof(values[MODELENDYEAR]);
1414 if (model->nMaxSecVar > 0)
1415 model->SecularVariationUsed = 1;
1416 else
1417 model->SecularVariationUsed = 0;
1418}
1419
1420void MAG_AssignMagneticModelCoeffs(MAGtype_MagneticModel *Assignee,
1421 MAGtype_MagneticModel *Source, int nMax,
1422 int nMaxSecVar)
1423/* This function assigns the first nMax degrees of the Source model to the
1424 Assignee model, leaving the other coefficients untouched*/
1425{
1426 int n, m, index;
1427 assert(nMax <= Source->nMax);
1428 assert(nMax <= Assignee->nMax);
1429 assert(nMaxSecVar <= Source->nMaxSecVar);
1430 assert(nMaxSecVar <= Assignee->nMaxSecVar);
1431 for (n = 1; n <= nMaxSecVar; n++) {
1432 for (m = 0; m <= n; m++) {
1433 index = (n * (n + 1) / 2 + m);
1434 Assignee->Main_Field_Coeff_G[index] = Source->Main_Field_Coeff_G[index];
1435 Assignee->Main_Field_Coeff_H[index] = Source->Main_Field_Coeff_H[index];
1436 Assignee->Secular_Var_Coeff_G[index] = Source->Secular_Var_Coeff_G[index];
1437 Assignee->Secular_Var_Coeff_H[index] = Source->Secular_Var_Coeff_H[index];
1438 }
1439 }
1440 for (n = nMaxSecVar + 1; n <= nMax; n++) {
1441 for (m = 0; m <= n; m++) {
1442 index = (n * (n + 1) / 2 + m);
1443 Assignee->Main_Field_Coeff_G[index] = Source->Main_Field_Coeff_G[index];
1444 Assignee->Main_Field_Coeff_H[index] = Source->Main_Field_Coeff_H[index];
1445 }
1446 }
1447 return;
1448} /*MAG_AssignMagneticModelCoeffs*/
1449
1450int MAG_FreeMemory(MAGtype_MagneticModel *MagneticModel,
1451 MAGtype_MagneticModel *TimedMagneticModel,
1452 MAGtype_LegendreFunction *LegendreFunction)
1453
1454/* Free memory used by WMM functions. Only to be called at the end of the main
1455function. INPUT : MagneticModel pointer to data structure with the following
1456elements
1457
1458 double EditionDate;
1459 double epoch; Base time of Geomagnetic model epoch
1460(yrs) char ModelName[20]; double *Main_Field_Coeff_G; C - Gauss
1461coefficients of main geomagnetic model (nT) double *Main_Field_Coeff_H; C -
1462Gauss coefficients of main geomagnetic model (nT) double *Secular_Var_Coeff_G;
1463CD - Gauss coefficients of secular geomagnetic model (nT/yr) double
1464*Secular_Var_Coeff_H; CD - Gauss coefficients of secular geomagnetic model
1465(nT/yr) int nMax; Maximum degree of spherical harmonic model int nMaxSecVar;
1466Maxumum degree of spherical harmonic secular model int SecularVariationUsed;
1467Whether or not the magnetic secular variation vector will be needed by program
1468
1469 TimedMagneticModel Pointer to data structure similar to the
1470first input. LegendreFunction Pointer to data structure with the following
1471elements double *Pcup; ( pointer to store Legendre Function ) double *dPcup;
1472( pointer to store Derivative of Lagendre function )
1473
1474OUTPUT none
1475CALLS : none
1476
1477 */
1478{
1479 if (MagneticModel->Main_Field_Coeff_G) {
1480 free(MagneticModel->Main_Field_Coeff_G);
1481 MagneticModel->Main_Field_Coeff_G = NULL;
1482 }
1483 if (MagneticModel->Main_Field_Coeff_H) {
1484 free(MagneticModel->Main_Field_Coeff_H);
1485 MagneticModel->Main_Field_Coeff_H = NULL;
1486 }
1487 if (MagneticModel->Secular_Var_Coeff_G) {
1488 free(MagneticModel->Secular_Var_Coeff_G);
1489 MagneticModel->Secular_Var_Coeff_G = NULL;
1490 }
1491 if (MagneticModel->Secular_Var_Coeff_H) {
1492 free(MagneticModel->Secular_Var_Coeff_H);
1493 MagneticModel->Secular_Var_Coeff_H = NULL;
1494 }
1495 if (MagneticModel) {
1496 free(MagneticModel);
1497 MagneticModel = NULL;
1498 }
1499
1500 if (TimedMagneticModel->Main_Field_Coeff_G) {
1501 free(TimedMagneticModel->Main_Field_Coeff_G);
1502 TimedMagneticModel->Main_Field_Coeff_G = NULL;
1503 }
1504 if (TimedMagneticModel->Main_Field_Coeff_H) {
1505 free(TimedMagneticModel->Main_Field_Coeff_H);
1506 TimedMagneticModel->Main_Field_Coeff_H = NULL;
1507 }
1508 if (TimedMagneticModel->Secular_Var_Coeff_G) {
1509 free(TimedMagneticModel->Secular_Var_Coeff_G);
1510 TimedMagneticModel->Secular_Var_Coeff_G = NULL;
1511 }
1512 if (TimedMagneticModel->Secular_Var_Coeff_H) {
1513 free(TimedMagneticModel->Secular_Var_Coeff_H);
1514 TimedMagneticModel->Secular_Var_Coeff_H = NULL;
1515 }
1516
1517 if (TimedMagneticModel) {
1518 free(TimedMagneticModel);
1519 TimedMagneticModel = NULL;
1520 }
1521
1522 if (LegendreFunction->Pcup) {
1523 free(LegendreFunction->Pcup);
1524 LegendreFunction->Pcup = NULL;
1525 }
1526 if (LegendreFunction->dPcup) {
1527 free(LegendreFunction->dPcup);
1528 LegendreFunction->dPcup = NULL;
1529 }
1530 if (LegendreFunction) {
1531 free(LegendreFunction);
1532 LegendreFunction = NULL;
1533 }
1534
1535 return TRUE;
1536} /*MAG_FreeMemory */
1537
1538int MAG_FreeMagneticModelMemory(MAGtype_MagneticModel *MagneticModel)
1539
1540/* Free the magnetic model memory used by WMM functions.
1541INPUT : MagneticModel pointer to data structure with the following elements
1542
1543 double EditionDate;
1544 double epoch; Base time of Geomagnetic model epoch
1545(yrs) char ModelName[20]; double *Main_Field_Coeff_G; C - Gauss
1546coefficients of main geomagnetic model (nT) double *Main_Field_Coeff_H; C -
1547Gauss coefficients of main geomagnetic model (nT) double *Secular_Var_Coeff_G;
1548CD - Gauss coefficients of secular geomagnetic model (nT/yr) double
1549*Secular_Var_Coeff_H; CD - Gauss coefficients of secular geomagnetic model
1550(nT/yr) int nMax; Maximum degree of spherical harmonic model int nMaxSecVar;
1551Maxumum degree of spherical harmonic secular model int SecularVariationUsed;
1552Whether or not the magnetic secular variation vector will be needed by program
1553
1554OUTPUT none
1555CALLS : none
1556
1557 */
1558{
1559 if (MagneticModel->Main_Field_Coeff_G) {
1560 free(MagneticModel->Main_Field_Coeff_G);
1561 MagneticModel->Main_Field_Coeff_G = NULL;
1562 }
1563 if (MagneticModel->Main_Field_Coeff_H) {
1564 free(MagneticModel->Main_Field_Coeff_H);
1565 MagneticModel->Main_Field_Coeff_H = NULL;
1566 }
1567 if (MagneticModel->Secular_Var_Coeff_G) {
1568 free(MagneticModel->Secular_Var_Coeff_G);
1569 MagneticModel->Secular_Var_Coeff_G = NULL;
1570 }
1571 if (MagneticModel->Secular_Var_Coeff_H) {
1572 free(MagneticModel->Secular_Var_Coeff_H);
1573 MagneticModel->Secular_Var_Coeff_H = NULL;
1574 }
1575 if (MagneticModel) {
1576 free(MagneticModel);
1577 MagneticModel = NULL;
1578 }
1579
1580 return TRUE;
1581} /*MAG_FreeMagneticModelMemory */
1582
1583int MAG_FreeLegendreMemory(MAGtype_LegendreFunction *LegendreFunction)
1584
1585/* Free the Legendre Coefficients memory used by the WMM functions.
1586INPUT : LegendreFunction Pointer to data structure with the following elements
1587 double *Pcup; ( pointer to
1588store Legendre Function ) double *dPcup; ( pointer to store Derivative of
1589Lagendre function )
1590
1591OUTPUT: none
1592CALLS : none
1593
1594 */
1595{
1596 if (LegendreFunction->Pcup) {
1597 free(LegendreFunction->Pcup);
1598 LegendreFunction->Pcup = NULL;
1599 }
1600 if (LegendreFunction->dPcup) {
1601 free(LegendreFunction->dPcup);
1602 LegendreFunction->dPcup = NULL;
1603 }
1604 if (LegendreFunction) {
1605 free(LegendreFunction);
1606 LegendreFunction = NULL;
1607 }
1608
1609 return TRUE;
1610} /*MAG_FreeLegendreMemory */
1611
1612int MAG_FreeSphVarMemory(MAGtype_SphericalHarmonicVariables *SphVar)
1613
1614/* Free the Spherical Harmonic Variable memory used by the WMM functions.
1615INPUT : LegendreFunction Pointer to data structure with the following elements
1616 double *RelativeRadiusPower
1617 double *cos_mlambda
1618 double *sin_mlambda
1619 OUTPUT: none
1620 CALLS : none
1621 */
1622{
1623 if (SphVar->RelativeRadiusPower) {
1624 free(SphVar->RelativeRadiusPower);
1625 SphVar->RelativeRadiusPower = NULL;
1626 }
1627 if (SphVar->cos_mlambda) {
1628 free(SphVar->cos_mlambda);
1629 SphVar->cos_mlambda = NULL;
1630 }
1631 if (SphVar->sin_mlambda) {
1632 free(SphVar->sin_mlambda);
1633 SphVar->sin_mlambda = NULL;
1634 }
1635 if (SphVar) {
1636 free(SphVar);
1637 SphVar = NULL;
1638 }
1639
1640 return TRUE;
1641} /*MAG_FreeSphVarMemory*/
1642
1643#ifndef OPENCPN
1644void MAG_PrintWMMFormat(char *filename, MAGtype_MagneticModel *MagneticModel) {
1645 int index, n, m;
1646 FILE *OUT;
1647 MAGtype_Date Date;
1648 char Datestring[11];
1649
1650 Date.DecimalYear = MagneticModel->EditionDate;
1651 MAG_YearToDate(&Date);
1652 sprintf(Datestring, "%d/%d/%d", Date.Month, Date.Day, Date.Year);
1653 OUT = fopen(filename, "w");
1654 fprintf(OUT, " %.1f %s %s\n",
1655 MagneticModel->epoch, MagneticModel->ModelName, Datestring);
1656 for (n = 1; n <= MagneticModel->nMax; n++) {
1657 for (m = 0; m <= n; m++) {
1658 index = (n * (n + 1) / 2 + m);
1659 if (m != 0)
1660 fprintf(OUT, " %2d %2d %9.4f %9.4f %9.4f %9.4f\n", n, m,
1661 MagneticModel->Main_Field_Coeff_G[index],
1662 MagneticModel->Main_Field_Coeff_H[index],
1663 MagneticModel->Secular_Var_Coeff_G[index],
1664 MagneticModel->Secular_Var_Coeff_H[index]);
1665 else
1666 fprintf(OUT, " %2d %2d %9.4f %9.4f %9.4f %9.4f\n", n, m,
1667 MagneticModel->Main_Field_Coeff_G[index], 0.0,
1668 MagneticModel->Secular_Var_Coeff_G[index], 0.0);
1669 }
1670 }
1671 fclose(OUT);
1672} /*MAG_PrintWMMFormat*/
1673
1674#endif /* OPENCPN */
1675
1676#ifndef OPENCPN
1677void MAG_PrintEMMFormat(char *filename, char *filenameSV,
1678 MAGtype_MagneticModel *MagneticModel) {
1679 int index, n, m;
1680 FILE *OUT;
1681 MAGtype_Date Date;
1682 char Datestring[11];
1683
1684 Date.DecimalYear = MagneticModel->EditionDate;
1685 MAG_YearToDate(&Date);
1686 sprintf(Datestring, "%d/%d/%d", Date.Month, Date.Day, Date.Year);
1687 OUT = fopen(filename, "w");
1688 fprintf(OUT, " %.1f %s %s\n",
1689 MagneticModel->epoch, MagneticModel->ModelName, Datestring);
1690 for (n = 1; n <= MagneticModel->nMax; n++) {
1691 for (m = 0; m <= n; m++) {
1692 index = (n * (n + 1) / 2 + m);
1693 if (m != 0)
1694 fprintf(OUT, " %2d %2d %9.4f %9.4f\n", n, m,
1695 MagneticModel->Main_Field_Coeff_G[index],
1696 MagneticModel->Main_Field_Coeff_H[index]);
1697 else
1698 fprintf(OUT, " %2d %2d %9.4f %9.4f\n", n, m,
1699 MagneticModel->Main_Field_Coeff_G[index], 0.0);
1700 }
1701 }
1702 fclose(OUT);
1703 OUT = fopen(filenameSV, "w");
1704 for (n = 1; n <= MagneticModel->nMaxSecVar; n++) {
1705 for (m = 0; m <= n; m++) {
1706 index = (n * (n + 1) / 2 + m);
1707 if (m != 0)
1708 fprintf(OUT, " %2d %2d %9.4f %9.4f\n", n, m,
1709 MagneticModel->Secular_Var_Coeff_G[index],
1710 MagneticModel->Secular_Var_Coeff_H[index]);
1711 else
1712 fprintf(OUT, " %2d %2d %9.4f %9.4f\n", n, m,
1713 MagneticModel->Secular_Var_Coeff_G[index], 0.0);
1714 }
1715 }
1716 fclose(OUT);
1717 return;
1718} /*MAG_PrintEMMFormat*/
1719
1720#endif /* OPENCPN */
1721
1722#ifndef OPENCPN
1723void MAG_PrintSHDFFormat(char *filename,
1724 MAGtype_MagneticModel *(*MagneticModel)[]) {
1725 int epochs = 1;
1726 int i, n, m, index, epochRange;
1727 FILE *SHDF_file;
1728 SHDF_file = fopen(filename, "w");
1729 /*lines = (int)(UFM_DEGREE / 2.0 * (UFM_DEGREE + 3));*/
1730 for (i = 0; i < epochs; i++) {
1731 if (i < epochs - 1)
1732 epochRange = (*MagneticModel)[i + 1]->epoch - (*MagneticModel)[i]->epoch;
1733 else
1734 epochRange = (*MagneticModel)[i]->epoch - (*MagneticModel)[i - 1]->epoch;
1735 fprintf(SHDF_file,
1736 "%%SHDF 16695 Definitive Geomagnetic Reference Field Model "
1737 "Coefficient File\n");
1738 fprintf(SHDF_file, "%%ModelName: %s\n", (*MagneticModel)[i]->ModelName);
1739 fprintf(SHDF_file,
1740 "%%Publisher: International Association of Geomagnetism and "
1741 "Aeronomy (IAGA), Working Group V-Mod\n");
1742 fprintf(SHDF_file, "%%ReleaseDate: Some Number\n");
1743 fprintf(SHDF_file, "%%DataCutOFF: Some Other Number\n");
1744 fprintf(SHDF_file, "%%ModelStartYear: %d\n",
1745 (int)(*MagneticModel)[i]->epoch);
1746 fprintf(SHDF_file, "%%ModelEndYear: %d\n",
1747 (int)(*MagneticModel)[i]->epoch + epochRange);
1748 fprintf(SHDF_file, "%%Epoch: %.0f\n", (*MagneticModel)[i]->epoch);
1749 fprintf(SHDF_file, "%%IntStaticDeg: %d\n", (*MagneticModel)[i]->nMax);
1750 fprintf(SHDF_file, "%%IntSecVarDeg: %d\n", (*MagneticModel)[i]->nMaxSecVar);
1751 fprintf(SHDF_file, "%%ExtStaticDeg: 0\n");
1752 fprintf(SHDF_file, "%%ExtSecVarDeg: 0\n");
1753 fprintf(SHDF_file, "%%Normalization: Schmidt semi-normailized\n");
1754 fprintf(SHDF_file, "%%SpatBasFunc: spherical harmonics\n");
1755 fprintf(SHDF_file, "# To synthesize the field for a given date:\n");
1756 fprintf(SHDF_file,
1757 "# Use the sub-model of the epoch corresponding to each date\n");
1758 fprintf(SHDF_file,
1759 "#\n#\n#\n#\n# I/E, n, m, Gnm, Hnm, SV-Gnm, SV-Hnm\n#\n");
1760 n = 1;
1761 m = 0;
1762 for (n = 1; n <= (*MagneticModel)[i]->nMax; n++) {
1763 for (m = 0; m <= n; m++) {
1764 index = (n * (n + 1)) / 2 + m;
1765 if (i < epochs - 1) {
1766 if (m != 0)
1767 fprintf(SHDF_file, "I,%d,%d,%f,%f,%f,%f\n", n, m,
1768 (*MagneticModel)[i]->Main_Field_Coeff_G[index],
1769 (*MagneticModel)[i]->Main_Field_Coeff_H[index],
1770 (*MagneticModel)[i]->Secular_Var_Coeff_G[index],
1771 (*MagneticModel)[i]->Secular_Var_Coeff_H[index]);
1772 else
1773 fprintf(SHDF_file, "I,%d,%d,%f,,%f,\n", n, m,
1774 (*MagneticModel)[i]->Main_Field_Coeff_G[index],
1775 (*MagneticModel)[i]->Secular_Var_Coeff_G[index]);
1776 } else {
1777 if (m != 0)
1778 fprintf(SHDF_file, "I,%d,%d,%f,%f,%f,%f\n", n, m,
1779 (*MagneticModel)[i]->Main_Field_Coeff_G[index],
1780 (*MagneticModel)[i]->Main_Field_Coeff_H[index],
1781 (*MagneticModel)[i]->Secular_Var_Coeff_G[index],
1782 (*MagneticModel)[i]->Secular_Var_Coeff_H[index]);
1783 else
1784 fprintf(SHDF_file, "I,%d,%d,%f,,%f,\n", n, m,
1785 (*MagneticModel)[i]->Main_Field_Coeff_G[index],
1786 (*MagneticModel)[i]->Secular_Var_Coeff_G[index]);
1787 }
1788 }
1789 }
1790 }
1791} /*MAG_PrintSHDFFormat*/
1792
1793#endif /* OPENCPN */
1794
1795int MAG_readMagneticModel(char *filename,
1796 MAGtype_MagneticModel *MagneticModel) {
1797 /* READ WORLD Magnetic MODEL SPHERICAL HARMONIC COEFFICIENTS (WMM.cof)
1798 INPUT : filename
1799 MagneticModel : Pointer to the data structure with the following
1800 fields required as inputs nMax : Number of static coefficients UPDATES :
1801 MagneticModel : Pointer to the data structure with the following fields
1802 populated char *ModelName; double epoch; Base time of Geomagnetic
1803 model epoch (yrs) double *Main_Field_Coeff_G; C - Gauss
1804 coefficients of main geomagnetic model (nT) double *Main_Field_Coeff_H; C -
1805 Gauss coefficients of main geomagnetic model (nT) double
1806 *Secular_Var_Coeff_G; CD - Gauss coefficients of secular geomagnetic model
1807 (nT/yr) double *Secular_Var_Coeff_H; CD - Gauss coefficients of secular
1808 geomagnetic model (nT/yr) CALLS : none
1809
1810 */
1811
1812 FILE *MAG_COF_File;
1813 char c_str[81],
1814 c_new[5]; /*these strings are used to read a line from coefficient file*/
1815 int i, icomp, m, n, EOF_Flag = 0, index;
1816 double epoch, gnm, hnm, dgnm, dhnm;
1817 MAG_COF_File = fopen(filename, "r");
1818
1819 if (MAG_COF_File == NULL) {
1820 MAG_Error(20);
1821 return FALSE;
1822 /* should we have a standard error printing routine ?*/
1823 }
1824 MagneticModel->Main_Field_Coeff_H[0] = 0.0;
1825 MagneticModel->Main_Field_Coeff_G[0] = 0.0;
1826 MagneticModel->Secular_Var_Coeff_H[0] = 0.0;
1827 MagneticModel->Secular_Var_Coeff_G[0] = 0.0;
1828 fgets(c_str, 80, MAG_COF_File);
1829 sscanf(c_str, "%lf%s", &epoch, MagneticModel->ModelName);
1830 MagneticModel->epoch = epoch;
1831 while (EOF_Flag == 0) {
1832 if (NULL == fgets(c_str, 80, MAG_COF_File)) {
1833 break;
1834 }
1835 /* CHECK FOR LAST LINE IN FILE */
1836 for (i = 0; i < 4 && (c_str[i] != '\0'); i++) {
1837 c_new[i] = c_str[i];
1838 c_new[i + 1] = '\0';
1839 }
1840 icomp = strcmp("9999", c_new);
1841 if (icomp == 0) {
1842 EOF_Flag = 1;
1843 break;
1844 }
1845 /* END OF FILE NOT ENCOUNTERED, GET VALUES */
1846 sscanf(c_str, "%d%d%lf%lf%lf%lf", &n, &m, &gnm, &hnm, &dgnm, &dhnm);
1847 if (m <= n) {
1848 index = (n * (n + 1) / 2 + m);
1849 MagneticModel->Main_Field_Coeff_G[index] = gnm;
1850 MagneticModel->Secular_Var_Coeff_G[index] = dgnm;
1851 MagneticModel->Main_Field_Coeff_H[index] = hnm;
1852 MagneticModel->Secular_Var_Coeff_H[index] = dhnm;
1853 }
1854 }
1855
1856 fclose(MAG_COF_File);
1857 return TRUE;
1858} /*MAG_readMagneticModel*/
1859
1860int MAG_readMagneticModel_Large(char *filename, char *filenameSV,
1861 MAGtype_MagneticModel *MagneticModel)
1862
1863/* To read the high-degree model coefficients (for example, NGDC 720)
1864 INPUT : filename file name for static coefficients
1865 filenameSV file name for secular variation coefficients
1866
1867 MagneticModel : Pointer to the data structure with the
1868 following fields required as inputs nMaxSecVar : Number of secular variation
1869 coefficients nMax : Number of static coefficients UPDATES : MagneticModel :
1870 Pointer to the data structure with the following fields populated double
1871 epoch; Base time of Geomagnetic model epoch (yrs) double
1872 *Main_Field_Coeff_G; C - Gauss coefficients of main geomagnetic
1873 model (nT) double *Main_Field_Coeff_H; C - Gauss coefficients of
1874 main geomagnetic model (nT) double *Secular_Var_Coeff_G; CD - Gauss
1875 coefficients of secular geomagnetic model (nT/yr) double
1876 *Secular_Var_Coeff_H; CD - Gauss coefficients of secular geomagnetic model
1877 (nT/yr) CALLS : none
1878
1879 */
1880{
1881 FILE *MAG_COF_File;
1882 FILE *MAG_COFSV_File;
1883 char c_str[81], c_str2[81]; /* these strings are used to read a line from
1884 coefficient file */
1885 int i, m, n, index, a, b;
1886 double epoch, gnm, hnm, dgnm, dhnm;
1887 MAG_COF_File = fopen(filename, "r");
1888 MAG_COFSV_File = fopen(filenameSV, "r");
1889 if (MAG_COF_File == NULL || MAG_COFSV_File == NULL) {
1890 MAG_Error(20);
1891 return FALSE;
1892 }
1893 MagneticModel->Main_Field_Coeff_H[0] = 0.0;
1894 MagneticModel->Main_Field_Coeff_G[0] = 0.0;
1895 MagneticModel->Secular_Var_Coeff_H[0] = 0.0;
1896 MagneticModel->Secular_Var_Coeff_G[0] = 0.0;
1897 if (NULL == fgets(c_str, 80, MAG_COF_File)) {
1898 fclose(MAG_COF_File);
1899 fclose(MAG_COFSV_File);
1900 return FALSE;
1901 }
1902 sscanf(c_str, "%lf%s", &epoch, MagneticModel->ModelName);
1903 MagneticModel->epoch = epoch;
1904 a = CALCULATE_NUMTERMS(MagneticModel->nMaxSecVar);
1905 b = CALCULATE_NUMTERMS(MagneticModel->nMax);
1906 for (i = 0; i < a; i++) {
1907 if (NULL == fgets(c_str, 80, MAG_COF_File)) {
1908 fclose(MAG_COF_File);
1909 fclose(MAG_COFSV_File);
1910 return FALSE;
1911 }
1912 sscanf(c_str, "%d%d%lf%lf", &n, &m, &gnm, &hnm);
1913 if (NULL == fgets(c_str2, 80, MAG_COFSV_File)) {
1914 fclose(MAG_COF_File);
1915 fclose(MAG_COFSV_File);
1916 return FALSE;
1917 }
1918 sscanf(c_str2, "%d%d%lf%lf", &n, &m, &dgnm, &dhnm);
1919 if (m <= n) {
1920 index = (n * (n + 1) / 2 + m);
1921 MagneticModel->Main_Field_Coeff_G[index] = gnm;
1922 MagneticModel->Secular_Var_Coeff_G[index] = dgnm;
1923 MagneticModel->Main_Field_Coeff_H[index] = hnm;
1924 MagneticModel->Secular_Var_Coeff_H[index] = dhnm;
1925 }
1926 }
1927 for (i = a; i < b; i++) {
1928 if (NULL == fgets(c_str, 80, MAG_COF_File)) {
1929 fclose(MAG_COF_File);
1930 fclose(MAG_COFSV_File);
1931 return FALSE;
1932 }
1933 sscanf(c_str, "%d%d%lf%lf", &n, &m, &gnm, &hnm);
1934 if (m <= n) {
1935 index = (n * (n + 1) / 2 + m);
1936 MagneticModel->Main_Field_Coeff_G[index] = gnm;
1937 MagneticModel->Main_Field_Coeff_H[index] = hnm;
1938 }
1939 }
1940 if (MAG_COF_File != NULL && MAG_COFSV_File != NULL) {
1941 fclose(MAG_COF_File);
1942 fclose(MAG_COFSV_File);
1943 }
1944
1945 return TRUE;
1946} /*MAG_readMagneticModel_Large*/
1947
1948int MAG_readMagneticModel_SHDF(char *filename,
1949 MAGtype_MagneticModel *(*magneticmodels)[1])
1950/*
1951 * MAG_readMagneticModels - Read the Magnetic Models from an SHDF format file
1952 *
1953 * Input:
1954 * filename - Path to the SHDF format model file to be read
1955 * array_size - Max No of models to be read from the file
1956 *
1957 * Output:
1958 * magneticmodels[] - Array of magnetic models read from the file
1959 *
1960 * Return value:
1961 * Returns the number of models read from the file.
1962 * -2 implies that internal or external static degree was not found in the
1963 * file, hence memory cannot be allocated -1 implies some error during file
1964 * processing (I/O) 0 implies no models were read from the file if ReturnValue >
1965 * array_size then there were too many models in model file but only
1966 * <array_size> number were read . if ReturnValue <= array_size then the
1967 * function execution was successful.
1968 */
1969{
1970 char paramkeys[NOOFPARAMS][MAXLINELENGTH] = {
1971 "SHDF ", "ModelName: ", "Publisher: ", "ReleaseDate: ",
1972 "DataCutOff: ", "ModelStartYear: ", "ModelEndYear: ", "Epoch: ",
1973 "IntStaticDeg: ", "IntSecVarDeg: ", "ExtStaticDeg: ", "ExtSecVarDeg: ",
1974 "GeoMagRefRad: ", "Normalization: ", "SpatBasFunc: "};
1975
1976 char paramvalues[NOOFPARAMS][MAXLINELENGTH];
1977 char *line = (char *)malloc(MAXLINELENGTH);
1978 char *ptrreset;
1979 char paramvalue[MAXLINELENGTH];
1980 int paramvaluelength = 0;
1981 int paramkeylength = 0;
1982 int i = 0, j = 0;
1983 int newrecord = 1;
1984 int header_index = -1;
1985 int numterms;
1986 int tempint;
1987 int allocationflag = 0;
1988 char coefftype; /* Internal or External (I/E) */
1989 int array_size = 1;
1990
1991 /* For reading coefficients */
1992 int n, m;
1993 double gnm, hnm, dgnm, dhnm, cutoff;
1994 int index;
1995
1996 FILE *stream;
1997 ptrreset = line;
1998 stream = fopen(filename, READONLYMODE);
1999 if (stream == NULL) {
2000 perror("File open error");
2001 return header_index;
2002 }
2003
2004 /* Read records from the model file and store header information. */
2005 while (fgets(line, MAXLINELENGTH, stream) != NULL) {
2006 j++;
2007 if (strlen(MAG_Trim(line)) == 0) continue;
2008 if (*line == '%') {
2009 line++;
2010 if (newrecord) {
2011 if (header_index > -1) {
2012 MAG_AssignHeaderValues((*magneticmodels)[header_index], paramvalues);
2013 }
2014 header_index++;
2015 if (header_index >= array_size) {
2016 fprintf(
2017 stderr,
2018 "Header limit exceeded - too many models in model file. (%d)\n",
2019 header_index);
2020 return array_size + 1;
2021 }
2022 newrecord = 0;
2023 allocationflag = 0;
2024 }
2025 for (i = 0; i < NOOFPARAMS; i++) {
2026 paramkeylength = strlen(paramkeys[i]);
2027 if (!strncmp(line, paramkeys[i], paramkeylength)) {
2028 paramvaluelength = strlen(line) - paramkeylength;
2029 strncpy(paramvalue, line + paramkeylength, paramvaluelength);
2030 paramvalue[paramvaluelength] = '\0';
2031 strcpy(paramvalues[i], paramvalue);
2032 if (!strcmp(paramkeys[i], paramkeys[INTSTATICDEG]) ||
2033 !strcmp(paramkeys[i], paramkeys[EXTSTATICDEG])) {
2034 tempint = atoi(paramvalues[i]);
2035 if (tempint > 0 && allocationflag == 0) {
2036 numterms = CALCULATE_NUMTERMS(tempint);
2037 (*magneticmodels)[header_index] =
2038 MAG_AllocateModelMemory(numterms);
2039 /* model = (*magneticmodels)[header_index]; */
2040 allocationflag = 1;
2041 }
2042 }
2043 break;
2044 }
2045 }
2046 line--;
2047 } else if (*line == '#') {
2048 /* process comments */
2049
2050 } else if (sscanf(line, "%c,%d,%d", &coefftype, &n, &m) == 3) {
2051 if (m == 0) {
2052 sscanf(line, "%c,%d,%d,%lf,,%lf,", &coefftype, &n, &m, &gnm, &dgnm);
2053 hnm = 0;
2054 dhnm = 0;
2055 } else
2056 sscanf(line, "%c,%d,%d,%lf,%lf,%lf,%lf", &coefftype, &n, &m, &gnm, &hnm,
2057 &dgnm, &dhnm);
2058 newrecord = 1;
2059 if (!allocationflag) {
2060 fprintf(stderr,
2061 "Degree not found in model. Memory cannot be allocated.\n");
2062 return _DEGREE_NOT_FOUND;
2063 }
2064 if (m <= n) {
2065 index = (n * (n + 1) / 2 + m);
2066 (*magneticmodels)[header_index]->Main_Field_Coeff_G[index] = gnm;
2067 (*magneticmodels)[header_index]->Secular_Var_Coeff_G[index] = dgnm;
2068 (*magneticmodels)[header_index]->Main_Field_Coeff_H[index] = hnm;
2069 (*magneticmodels)[header_index]->Secular_Var_Coeff_H[index] = dhnm;
2070 }
2071 }
2072 }
2073 if (header_index > -1)
2074 MAG_AssignHeaderValues((*magneticmodels)[header_index], paramvalues);
2075 fclose(stream);
2076
2077 cutoff = (*magneticmodels)[array_size - 1]->CoefficientFileEndDate;
2078
2079 for (i = 0; i < array_size; i++)
2080 (*magneticmodels)[i]->CoefficientFileEndDate = cutoff;
2081
2082 free(ptrreset);
2083 line = NULL;
2084 ptrreset = NULL;
2085 return header_index + 1;
2086} /*MAG_readMagneticModel_SHDF*/
2087
2088char *MAG_Trim(char *str) {
2089 char *end;
2090
2091 while (isspace(*str)) str++;
2092
2093 if (*str == 0) return str;
2094
2095 end = str + strlen(str) - 1;
2096 while (end > str && isspace(*end)) end--;
2097
2098 *(end + 1) = 0;
2099
2100 return str;
2101}
2102
2103/*End of Memory and File Processing functions*/
2104
2105/******************************************************************************
2106 *************Conversions, Transformations, and other Calculations**************
2107 * This grouping consists of functions that perform unit conversions, coordinate
2108 * transformations and other simple or straightforward calculations that are
2109 * usually easily replicable with a typical scientific calculator.
2110 ******************************************************************************/
2111
2112void MAG_BaseErrors(double DeclCoef, double DeclBaseline, double InclOffset,
2113 double FOffset, double Multiplier, double H,
2114 double *DeclErr, double *InclErr, double *FErr) {
2115 double declHorizontalAdjustmentSq;
2116 declHorizontalAdjustmentSq = (DeclCoef / H) * (DeclCoef / H);
2117 *DeclErr = sqrt(declHorizontalAdjustmentSq + DeclBaseline * DeclBaseline) *
2118 Multiplier;
2119 *InclErr = InclOffset * Multiplier;
2120 *FErr = FOffset * Multiplier;
2121}
2122
2123int MAG_CalculateGeoMagneticElements(
2124 MAGtype_MagneticResults *MagneticResultsGeo,
2125 MAGtype_GeoMagneticElements *GeoMagneticElements)
2126
2127/* Calculate all the Geomagnetic elements from X,Y and Z components
2128INPUT MagneticResultsGeo Pointer to data structure with the following
2129elements double Bx; ( North ) double By; ( East ) double Bz; ( Down )
2130OUTPUT GeoMagneticElements Pointer to data structure with the following
2131elements double Decl; (Angle between the magnetic field vector and true north,
2132positive east) double Incl; Angle between the magnetic field vector and the
2133horizontal plane, positive down double F; Magnetic Field Strength double H;
2134Horizontal Magnetic Field Strength double X; Northern component of the magnetic
2135field vector double Y; Eastern component of the magnetic field vector double Z;
2136Downward component of the magnetic field vector CALLS : none
2137 */
2138{
2139 GeoMagneticElements->X = MagneticResultsGeo->Bx;
2140 GeoMagneticElements->Y = MagneticResultsGeo->By;
2141 GeoMagneticElements->Z = MagneticResultsGeo->Bz;
2142
2143 GeoMagneticElements->H =
2144 sqrt(MagneticResultsGeo->Bx * MagneticResultsGeo->Bx +
2145 MagneticResultsGeo->By * MagneticResultsGeo->By);
2146 GeoMagneticElements->F =
2147 sqrt(GeoMagneticElements->H * GeoMagneticElements->H +
2148 MagneticResultsGeo->Bz * MagneticResultsGeo->Bz);
2149 GeoMagneticElements->Decl =
2150 RAD2DEG(atan2(GeoMagneticElements->Y, GeoMagneticElements->X));
2151 GeoMagneticElements->Incl =
2152 RAD2DEG(atan2(GeoMagneticElements->Z, GeoMagneticElements->H));
2153
2154 return TRUE;
2155} /*MAG_CalculateGeoMagneticElements */
2156
2157int MAG_CalculateGridVariation(MAGtype_CoordGeodetic location,
2158 MAGtype_GeoMagneticElements *elements)
2159
2160/*Computes the grid variation for |latitudes| > MAG_MAX_LAT_DEGREE
2161
2162Grivation (or grid variation) is the angle between grid north and
2163magnetic north. This routine calculates Grivation for the Polar Stereographic
2164projection for polar locations (Latitude => |55| deg). Otherwise, it computes
2165the grid variation in UTM projection system. However, the UTM projection codes
2166may be used to compute the grid variation at any latitudes.
2167
2168INPUT location Data structure with the following elements
2169 double lambda; (longitude)
2170 double phi; ( geodetic latitude)
2171 double HeightAboveEllipsoid; (height above the ellipsoid (HaE) )
2172 double HeightAboveGeoid;(height above the Geoid )
2173OUTPUT elements Data structure with the following elements updated
2174 double GV; ( The Grid Variation )
2175CALLS : MAG_GetTransverseMercator
2176
2177 */
2178{
2179 MAGtype_UTMParameters UTMParameters;
2180 if (location.phi >= MAG_PS_MAX_LAT_DEGREE) {
2181 elements->GV = elements->Decl - location.lambda;
2182 return 1;
2183 } else if (location.phi <= MAG_PS_MIN_LAT_DEGREE) {
2184 elements->GV = elements->Decl + location.lambda;
2185 return 1;
2186 } else {
2187 MAG_GetTransverseMercator(location, &UTMParameters);
2188 elements->GV = elements->Decl - UTMParameters.ConvergenceOfMeridians;
2189 }
2190 return 0;
2191} /*MAG_CalculateGridVariation*/
2192
2193void MAG_CalculateGradientElements(MAGtype_MagneticResults GradResults,
2194 MAGtype_GeoMagneticElements MagneticElements,
2195 MAGtype_GeoMagneticElements *GradElements) {
2196 GradElements->X = GradResults.Bx;
2197 GradElements->Y = GradResults.By;
2198 GradElements->Z = GradResults.Bz;
2199
2200 GradElements->H = (GradElements->X * MagneticElements.X +
2201 GradElements->Y * MagneticElements.Y) /
2202 MagneticElements.H;
2203 GradElements->F = (GradElements->X * MagneticElements.X +
2204 GradElements->Y * MagneticElements.Y +
2205 GradElements->Z * MagneticElements.Z) /
2206 MagneticElements.F;
2207 GradElements->Decl = 180.0 / M_PI *
2208 (MagneticElements.X * GradElements->Y -
2209 MagneticElements.Y * GradElements->X) /
2210 (MagneticElements.H * MagneticElements.H);
2211 GradElements->Incl = 180.0 / M_PI *
2212 (MagneticElements.H * GradElements->Z -
2213 MagneticElements.Z * GradElements->H) /
2214 (MagneticElements.F * MagneticElements.F);
2215 GradElements->GV = GradElements->Decl;
2216}
2217
2218int MAG_CalculateSecularVariationElements(
2219 MAGtype_MagneticResults MagneticVariation,
2220 MAGtype_GeoMagneticElements *MagneticElements)
2221/*This takes the Magnetic Variation in x, y, and z and uses it to calculate the
2222 secular variation of each of the Geomagnetic elements. INPUT
2223 MagneticVariation Data structure with the following elements double Bx; (
2224 North ) double By; ( East ) double Bz; ( Down ) OUTPUT MagneticElements
2225 Pointer to the data structure with the following elements updated double
2226 Decldot; Yearly Rate of change in declination double Incldot; Yearly Rate of
2227 change in inclination double Fdot; Yearly rate of change in Magnetic field
2228 strength double Hdot; Yearly rate of change in horizontal field strength
2229 double Xdot; Yearly rate of change in the northern
2230 component double Ydot; Yearly rate of change in the eastern component double
2231 Zdot; Yearly rate of change in the downward component double GVdot;Yearly
2232 rate of chnage in grid variation CALLS : none
2233
2234 */
2235{
2236 MagneticElements->Xdot = MagneticVariation.Bx;
2237 MagneticElements->Ydot = MagneticVariation.By;
2238 MagneticElements->Zdot = MagneticVariation.Bz;
2239 MagneticElements->Hdot =
2240 (MagneticElements->X * MagneticElements->Xdot +
2241 MagneticElements->Y * MagneticElements->Ydot) /
2242 MagneticElements->H; /* See equation 19 in the WMM technical report */
2243 MagneticElements->Fdot = (MagneticElements->X * MagneticElements->Xdot +
2244 MagneticElements->Y * MagneticElements->Ydot +
2245 MagneticElements->Z * MagneticElements->Zdot) /
2246 MagneticElements->F;
2247 MagneticElements->Decldot = 180.0 / M_PI *
2248 (MagneticElements->X * MagneticElements->Ydot -
2249 MagneticElements->Y * MagneticElements->Xdot) /
2250 (MagneticElements->H * MagneticElements->H);
2251 MagneticElements->Incldot = 180.0 / M_PI *
2252 (MagneticElements->H * MagneticElements->Zdot -
2253 MagneticElements->Z * MagneticElements->Hdot) /
2254 (MagneticElements->F * MagneticElements->F);
2255 MagneticElements->GVdot = MagneticElements->Decldot;
2256 return TRUE;
2257} /*MAG_CalculateSecularVariationElements*/
2258
2259void MAG_CartesianToGeodetic(MAGtype_Ellipsoid Ellip, double x, double y,
2260 double z, MAGtype_CoordGeodetic *CoordGeodetic) {
2261 /*This converts the Cartesian x, y, and z coordinates to Geodetic Coordinates
2262 x is defined as the direction pointing out of the core toward the point
2263 defined
2264 * by 0 degrees latitude and longitude.
2265 y is defined as the direction from the core toward 90 degrees east longitude
2266 along
2267 * the equator
2268 z is defined as the direction from the core out the geographic north pole*/
2269
2270 double modified_b, r, e, f, p, q, d, v, g, t, zlong, rlat;
2271
2272 /*
2273 * 1.0 compute semi-minor axis and set sign to that of z in order
2274 * to get sign of Phi correct
2275 */
2276
2277 if (z < 0.0)
2278 modified_b = -Ellip.b;
2279 else
2280 modified_b = Ellip.b;
2281
2282 /*
2283 * 2.0 compute intermediate values for latitude
2284 */
2285 r = sqrt(x * x + y * y);
2286 e = (modified_b * z - (Ellip.a * Ellip.a - modified_b * modified_b)) /
2287 (Ellip.a * r);
2288 f = (modified_b * z + (Ellip.a * Ellip.a - modified_b * modified_b)) /
2289 (Ellip.a * r);
2290 /*
2291 * 3.0 find solution to:
2292 * t^4 + 2*E*t^3 + 2*F*t - 1 = 0
2293 */
2294 p = (4.0 / 3.0) * (e * f + 1.0);
2295 q = 2.0 * (e * e - f * f);
2296 d = p * p * p + q * q;
2297
2298 if (d >= 0.0) {
2299 v = pow((sqrt(d) - q), (1.0 / 3.0)) - pow((sqrt(d) + q), (1.0 / 3.0));
2300 } else {
2301 v = 2.0 * sqrt(-p) * cos(acos(q / (p * sqrt(-p))) / 3.0);
2302 }
2303 /*
2304 * 4.0 improve v
2305 * NOTE: not really necessary unless point is near pole
2306 */
2307 if (v * v < fabs(p)) {
2308 v = -(v * v * v + 2.0 * q) / (3.0 * p);
2309 }
2310 g = (sqrt(e * e + v) + e) / 2.0;
2311 t = sqrt(g * g + (f - v * g) / (2.0 * g - e)) - g;
2312
2313 rlat = atan((Ellip.a * (1.0 - t * t)) / (2.0 * modified_b * t));
2314 CoordGeodetic->phi = RAD2DEG(rlat);
2315
2316 /*
2317 * 5.0 compute height above ellipsoid
2318 */
2319 CoordGeodetic->HeightAboveEllipsoid =
2320 (r - Ellip.a * t) * cos(rlat) + (z - modified_b) * sin(rlat);
2321 /*
2322 * 6.0 compute longitude east of Greenwich
2323 */
2324 zlong = atan2(y, x);
2325 if (zlong < 0.0) zlong = zlong + 2 * M_PI;
2326
2327 CoordGeodetic->lambda = RAD2DEG(zlong);
2328 while (CoordGeodetic->lambda > 180) {
2329 CoordGeodetic->lambda -= 360;
2330 }
2331}
2332
2333MAGtype_CoordGeodetic MAG_CoordGeodeticAssign(
2334 MAGtype_CoordGeodetic CoordGeodetic) {
2335 MAGtype_CoordGeodetic Assignee;
2336 Assignee.phi = CoordGeodetic.phi;
2337 Assignee.lambda = CoordGeodetic.lambda;
2338 Assignee.HeightAboveEllipsoid = CoordGeodetic.HeightAboveEllipsoid;
2339 Assignee.HeightAboveGeoid = CoordGeodetic.HeightAboveGeoid;
2340 Assignee.UseGeoid = CoordGeodetic.UseGeoid;
2341 return Assignee;
2342}
2343
2344int MAG_DateToYear(MAGtype_Date *CalendarDate, char *Error)
2345
2346/* Converts a given calendar date into a decimal year,
2347it also outputs an error string if there is a problem
2348INPUT CalendarDate Pointer to the data structure with the following elements
2349 int Year;
2350 int Month;
2351 int Day;
2352 double DecimalYear; decimal years
2353OUTPUT CalendarDate Pointer to the data structure with the following
2354elements updated double DecimalYear; decimal years Error pointer to an
2355error string CALLS : none
2356
2357 */
2358{
2359 int temp = 0; /*Total number of days */
2360 int MonthDays[13];
2361 int ExtraDay = 0;
2362 int i;
2363 if (CalendarDate->Month == 0) {
2364 CalendarDate->DecimalYear = CalendarDate->Year;
2365 return TRUE;
2366 }
2367 if ((CalendarDate->Year % 4 == 0 && CalendarDate->Year % 100 != 0) ||
2368 CalendarDate->Year % 400 == 0)
2369 ExtraDay = 1;
2370 MonthDays[0] = 0;
2371 MonthDays[1] = 31;
2372 MonthDays[2] = 28 + ExtraDay;
2373 MonthDays[3] = 31;
2374 MonthDays[4] = 30;
2375 MonthDays[5] = 31;
2376 MonthDays[6] = 30;
2377 MonthDays[7] = 31;
2378 MonthDays[8] = 31;
2379 MonthDays[9] = 30;
2380 MonthDays[10] = 31;
2381 MonthDays[11] = 30;
2382 MonthDays[12] = 31;
2383
2384 /******************Validation********************************/
2385 if (CalendarDate->Month <= 0 || CalendarDate->Month > 12) {
2386 strcpy(
2387 Error,
2388 "\nError: The Month entered is invalid, valid months are '1 to 12'\n");
2389 return 0;
2390 }
2391 if (CalendarDate->Day <= 0 ||
2392 CalendarDate->Day > MonthDays[CalendarDate->Month]) {
2393 printf("\nThe number of days in month %d is %d\n", CalendarDate->Month,
2394 MonthDays[CalendarDate->Month]);
2395 strcpy(Error, "\nError: The day entered is invalid\n");
2396 return 0;
2397 }
2398 /****************Calculation of t***************************/
2399 for (i = 1; i <= CalendarDate->Month; i++) temp += MonthDays[i - 1];
2400 temp += CalendarDate->Day;
2401 CalendarDate->DecimalYear =
2402 CalendarDate->Year + (temp - 1) / (365.0 + ExtraDay);
2403 return TRUE;
2404
2405} /*MAG_DateToYear*/
2406
2407void MAG_DegreeToDMSstring(double DegreesOfArc, int UnitDepth, char *DMSstring)
2408
2409/*This converts a given decimal degree into a DMS string.
2410INPUT DegreesOfArc decimal degree
2411 UnitDepth How many iterations should be printed,
2412 1 = Degrees
2413 2 = Degrees, Minutes
2414 3 = Degrees, Minutes, Seconds
2415OUPUT DMSstring pointer to DMSString. Must be at least 30 characters.
2416CALLS : none
2417 */
2418{
2419 int DMS[3], i;
2420 double temp = DegreesOfArc;
2421 char tempstring[36] = "";
2422 char tempstring2[32] = "";
2423 strcpy(DMSstring, "");
2424 if (UnitDepth > 3) MAG_Error(21);
2425 for (i = 0; i < UnitDepth; i++) {
2426 DMS[i] = (int)temp;
2427 switch (i) {
2428 case 0:
2429 strcpy(tempstring2, "Deg");
2430 break;
2431 case 1:
2432 strcpy(tempstring2, "Min");
2433 break;
2434 case 2:
2435 strcpy(tempstring2, "Sec");
2436 break;
2437 }
2438 temp = (temp - DMS[i]) * 60;
2439 if (i == UnitDepth - 1 && temp >= 30)
2440 DMS[i]++;
2441 else if (i == UnitDepth - 1 && temp <= -30)
2442 DMS[i]--;
2443 sprintf(tempstring, "%4d%4s", DMS[i], tempstring2);
2444 strcat(DMSstring, tempstring);
2445 }
2446} /*MAG_DegreeToDMSstring*/
2447
2448void MAG_DMSstringToDegree(char *DMSstring, double *DegreesOfArc)
2449
2450/*This converts a given DMS string into decimal degrees.
2451INPUT DMSstring pointer to DMSString
2452OUTPUT DegreesOfArc decimal degree
2453CALLS : none
2454 */
2455{
2456 int second, minute, degree, sign = 1, j = 0;
2457 j = sscanf(DMSstring, "%d, %d, %d", &degree, &minute, &second);
2458 if (j != 3) sscanf(DMSstring, "%d %d %d", &degree, &minute, &second);
2459 if (degree < 0) sign = -1;
2460 degree = degree * sign;
2461 *DegreesOfArc = sign * (degree + minute / 60.0 + second / 3600.0);
2462} /*MAG_DMSstringToDegree*/
2463
2464void MAG_ErrorCalc(MAGtype_GeoMagneticElements B,
2465 MAGtype_GeoMagneticElements *Errors) {
2466 /*Errors.Decl, Errors.Incl, Errors.F are all assumed to exist*/
2467 double cos2D, cos2I, sin2D, sin2I, EDSq, EISq, eD, eI;
2468 cos2D = cos(DEG2RAD(B.Decl)) * cos(DEG2RAD(B.Decl));
2469 cos2I = cos(DEG2RAD(B.Incl)) * cos(DEG2RAD(B.Incl));
2470 sin2D = sin(DEG2RAD(B.Decl)) * sin(DEG2RAD(B.Decl));
2471 sin2I = sin(DEG2RAD(B.Incl)) * sin(DEG2RAD(B.Incl));
2472 eD = DEG2RAD(Errors->Decl);
2473 eI = DEG2RAD(Errors->Incl);
2474 EDSq = eD * eD;
2475 EISq = eI * eI;
2476 Errors->X =
2477 sqrt(cos2D * cos2I * Errors->F * Errors->F +
2478 B.F * B.F * sin2D * cos2I * EDSq + B.F * B.F * cos2D * sin2I * EISq);
2479 Errors->Y =
2480 sqrt(sin2D * cos2I * Errors->F * Errors->F +
2481 B.F * B.F * cos2D * cos2I * EDSq + B.F * B.F * sin2D * sin2I * EISq);
2482 Errors->Z = sqrt(sin2I * Errors->F * Errors->F + B.F * B.F * cos2I * EISq);
2483 Errors->H = sqrt(cos2I * Errors->F * Errors->F + B.F * B.F * sin2I * EISq);
2484}
2485
2486int MAG_GeodeticToSpherical(MAGtype_Ellipsoid Ellip,
2487 MAGtype_CoordGeodetic CoordGeodetic,
2488 MAGtype_CoordSpherical *CoordSpherical)
2489
2490/* Converts Geodetic coordinates to Spherical coordinates
2491
2492 INPUT Ellip data structure with the following elements
2493 double a; semi-major axis of the ellipsoid
2494 double b; semi-minor axis of the ellipsoid
2495 double fla; flattening
2496 double epssq; first eccentricity squared
2497 double eps; first eccentricity
2498 double re; mean radius of ellipsoid
2499
2500 CoordGeodetic Pointer to the data structure with the
2501following elements updates double lambda; ( longitude ) double phi; ( geodetic
2502latitude ) double HeightAboveEllipsoid; ( height above the WGS84 ellipsoid (HaE)
2503) double HeightAboveGeoid; (height above the EGM96 Geoid model )
2504
2505 OUTPUT CoordSpherical Pointer to the data structure with the following
2506elements double lambda; ( longitude) double phig; ( geocentric latitude ) double
2507r; ( distance from the center of the ellipsoid)
2508
2509CALLS : none
2510
2511 */
2512{
2513 double CosLat, SinLat, rc, xp, zp; /*all local variables */
2514
2515 /*
2516 ** Convert geodetic coordinates, (defined by the WGS-84
2517 ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
2518 ** coordinates, and then to spherical coordinates.
2519 */
2520
2521 CosLat = cos(DEG2RAD(CoordGeodetic.phi));
2522 SinLat = sin(DEG2RAD(CoordGeodetic.phi));
2523
2524 /* compute the local radius of curvature on the WGS-84 reference ellipsoid */
2525
2526 rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);
2527
2528 /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */
2529
2530 xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;
2531 zp = (rc * (1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;
2532
2533 /* compute spherical radius and angle lambda and phi of specified point */
2534
2535 CoordSpherical->r = sqrt(xp * xp + zp * zp);
2536 CoordSpherical->phig =
2537 RAD2DEG(asin(zp / CoordSpherical->r)); /* geocentric latitude */
2538 CoordSpherical->lambda = CoordGeodetic.lambda; /* longitude */
2539
2540 return TRUE;
2541} /*MAG_GeodeticToSpherical*/
2542
2543MAGtype_GeoMagneticElements MAG_GeoMagneticElementsAssign(
2544 MAGtype_GeoMagneticElements Elements) {
2545 MAGtype_GeoMagneticElements Assignee;
2546 Assignee.X = Elements.X;
2547 Assignee.Y = Elements.Y;
2548 Assignee.Z = Elements.Z;
2549 Assignee.H = Elements.H;
2550 Assignee.F = Elements.F;
2551 Assignee.Decl = Elements.Decl;
2552 Assignee.Incl = Elements.Incl;
2553 Assignee.GV = Elements.GV;
2554 Assignee.Xdot = Elements.Xdot;
2555 Assignee.Ydot = Elements.Ydot;
2556 Assignee.Zdot = Elements.Zdot;
2557 Assignee.Hdot = Elements.Hdot;
2558 Assignee.Fdot = Elements.Fdot;
2559 Assignee.Decldot = Elements.Decldot;
2560 Assignee.Incldot = Elements.Incldot;
2561 Assignee.GVdot = Elements.GVdot;
2562 return Assignee;
2563}
2564
2565MAGtype_GeoMagneticElements MAG_GeoMagneticElementsScale(
2566 MAGtype_GeoMagneticElements Elements, double factor) {
2567 /*This function scales all the geomagnetic elements to scale a vector use
2568 MAG_MagneticResultsScale*/
2569 MAGtype_GeoMagneticElements product;
2570 product.X = Elements.X * factor;
2571 product.Y = Elements.Y * factor;
2572 product.Z = Elements.Z * factor;
2573 product.H = Elements.H * factor;
2574 product.F = Elements.F * factor;
2575 product.Incl = Elements.Incl * factor;
2576 product.Decl = Elements.Decl * factor;
2577 product.GV = Elements.GV * factor;
2578 product.Xdot = Elements.Xdot * factor;
2579 product.Ydot = Elements.Ydot * factor;
2580 product.Zdot = Elements.Zdot * factor;
2581 product.Hdot = Elements.Hdot * factor;
2582 product.Fdot = Elements.Fdot * factor;
2583 product.Incldot = Elements.Incldot * factor;
2584 product.Decldot = Elements.Decldot * factor;
2585 product.GVdot = Elements.GVdot * factor;
2586 return product;
2587}
2588
2589MAGtype_GeoMagneticElements MAG_GeoMagneticElementsSubtract(
2590 MAGtype_GeoMagneticElements minuend,
2591 MAGtype_GeoMagneticElements subtrahend) {
2592 /*This algorithm does not result in the difference of F being derived from
2593 the Pythagorean theorem. This function should be used for computing
2594 residuals or changes in elements.*/
2595 MAGtype_GeoMagneticElements difference;
2596 difference.X = minuend.X - subtrahend.X;
2597 difference.Y = minuend.Y - subtrahend.Y;
2598 difference.Z = minuend.Z - subtrahend.Z;
2599
2600 difference.H = minuend.H - subtrahend.H;
2601 difference.F = minuend.F - subtrahend.F;
2602 difference.Decl = minuend.Decl - subtrahend.Decl;
2603 difference.Incl = minuend.Incl - subtrahend.Incl;
2604
2605 difference.Xdot = minuend.Xdot - subtrahend.Xdot;
2606 difference.Ydot = minuend.Ydot - subtrahend.Ydot;
2607 difference.Zdot = minuend.Zdot - subtrahend.Zdot;
2608
2609 difference.Hdot = minuend.Hdot - subtrahend.Hdot;
2610 difference.Fdot = minuend.Fdot - subtrahend.Fdot;
2611 difference.Decldot = minuend.Decldot - subtrahend.Decldot;
2612 difference.Incldot = minuend.Incldot - subtrahend.Incldot;
2613
2614 difference.GV = minuend.GV - subtrahend.GV;
2615 difference.GVdot = minuend.GVdot - subtrahend.GVdot;
2616
2617 return difference;
2618}
2619
2620int MAG_GetTransverseMercator(MAGtype_CoordGeodetic CoordGeodetic,
2621 MAGtype_UTMParameters *UTMParameters)
2622/* Gets the UTM Parameters for a given Latitude and Longitude.
2623
2624INPUT: CoordGeodetic : Data structure MAGtype_CoordGeodetic.
2625OUTPUT : UTMParameters : Pointer to data structure MAGtype_UTMParameters with
2626the following elements double Easting; (X) in meters double Northing; (Y) in
2627meters int Zone; UTM Zone char HemiSphere ; double CentralMeridian; Longitude of
2628the Central Meridian of the UTM Zone double ConvergenceOfMeridians; Convergence
2629of Meridians double PointScale;
2630 */
2631{
2632 double Eps, Epssq;
2633 double Acoeff[8];
2634 double Lam0, K0, falseE, falseN;
2635 double K0R4, K0R4oa;
2636 double Lambda, Phi;
2637 int XYonly;
2638 double X, Y, pscale, CoM;
2639 int Zone;
2640 char Hemisphere;
2641
2642 /* Get the map projection parameters */
2643
2644 Lambda = DEG2RAD(CoordGeodetic.lambda);
2645 Phi = DEG2RAD(CoordGeodetic.phi);
2646
2647 MAG_GetUTMParameters(Phi, Lambda, &Zone, &Hemisphere, &Lam0);
2648 K0 = 0.9996;
2649
2650 if (Hemisphere == 'n' || Hemisphere == 'N') {
2651 falseN = 0;
2652 }
2653 if (Hemisphere == 's' || Hemisphere == 'S') {
2654 falseN = 10000000;
2655 }
2656
2657 falseE = 500000;
2658
2659 /* WGS84 ellipsoid */
2660
2661 Eps = 0.081819190842621494335;
2662 Epssq = 0.0066943799901413169961;
2663 K0R4 = 6367449.1458234153093 * K0;
2664 K0R4oa = K0R4 / 6378137;
2665
2666 Acoeff[0] = 8.37731820624469723600E-04;
2667 Acoeff[1] = 7.60852777357248641400E-07;
2668 Acoeff[2] = 1.19764550324249124400E-09;
2669 Acoeff[3] = 2.42917068039708917100E-12;
2670 Acoeff[4] = 5.71181837042801392800E-15;
2671 Acoeff[5] = 1.47999793137966169400E-17;
2672 Acoeff[6] = 4.10762410937071532000E-20;
2673 Acoeff[7] = 1.21078503892257704200E-22;
2674
2675 /* WGS84 ellipsoid */
2676
2677 /* Execution of the forward T.M. algorithm */
2678
2679 XYonly = 0;
2680
2681 MAG_TMfwd4(Eps, Epssq, K0R4, K0R4oa, Acoeff, Lam0, K0, falseE, falseN, XYonly,
2682 Lambda, Phi, &X, &Y, &pscale, &CoM);
2683
2684 /* Report results */
2685
2686 UTMParameters->Easting = X; /* UTM Easting (X) in meters*/
2687 UTMParameters->Northing = Y; /* UTM Northing (Y) in meters */
2688 UTMParameters->Zone = Zone; /*UTM Zone*/
2689 UTMParameters->HemiSphere = Hemisphere;
2690 UTMParameters->CentralMeridian =
2691 RAD2DEG(Lam0); /* Central Meridian of the UTM Zone */
2692 UTMParameters->ConvergenceOfMeridians =
2693 RAD2DEG(CoM); /* Convergence of meridians of the UTM Zone and location */
2694 UTMParameters->PointScale = pscale;
2695
2696 return 0;
2697} /*MAG_GetTransverseMercator*/
2698
2699int MAG_GetUTMParameters(double Latitude, double Longitude, int *Zone,
2700 char *Hemisphere, double *CentralMeridian) {
2701 /*
2702 * The function MAG_GetUTMParameters converts geodetic (latitude and
2703 * longitude) coordinates to UTM projection parameters (zone, hemisphere and
2704 * central meridian) If any errors occur, the error code(s) are returned by
2705 * the function, otherwise TRUE is returned.
2706 *
2707 * Latitude : Latitude in radians (input)
2708 * Longitude : Longitude in radians (input)
2709 * Zone : UTM zone (output)
2710 * Hemisphere : North or South hemisphere (output)
2711 * CentralMeridian : Central Meridian of the UTM Zone in radians (output)
2712 */
2713
2714 long Lat_Degrees;
2715 long Long_Degrees;
2716 long temp_zone;
2717 int Error_Code = 0;
2718
2719 if ((Latitude < DEG2RAD(MAG_UTM_MIN_LAT_DEGREE)) ||
2720 (Latitude >
2721 DEG2RAD(MAG_UTM_MAX_LAT_DEGREE))) { /* Latitude out of range */
2722 MAG_Error(23);
2723 Error_Code = 1;
2724 }
2725 if ((Longitude < -M_PI) ||
2726 (Longitude > (2 * M_PI))) { /* Longitude out of range */
2727 MAG_Error(24);
2728 Error_Code = 1;
2729 }
2730 if (!Error_Code) { /* no errors */
2731 if (Longitude < 0) Longitude += (2 * M_PI) + 1.0e-10;
2732 Lat_Degrees = (long)(Latitude * 180.0 / M_PI);
2733 Long_Degrees = (long)(Longitude * 180.0 / M_PI);
2734
2735 if (Longitude < M_PI)
2736 temp_zone = (long)(31 + ((Longitude * 180.0 / M_PI) / 6.0));
2737 else
2738 temp_zone = (long)(((Longitude * 180.0 / M_PI) / 6.0) - 29);
2739 if (temp_zone > 60) temp_zone = 1;
2740 /* UTM special cases */
2741 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1) &&
2742 (Long_Degrees < 3))
2743 temp_zone = 31;
2744 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2) &&
2745 (Long_Degrees < 12))
2746 temp_zone = 32;
2747 if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
2748 temp_zone = 31;
2749 if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
2750 temp_zone = 33;
2751 if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
2752 temp_zone = 35;
2753 if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
2754 temp_zone = 37;
2755
2756 if (!Error_Code) {
2757 if (temp_zone >= 31)
2758 *CentralMeridian = (6 * temp_zone - 183) * M_PI / 180.0;
2759 else
2760 *CentralMeridian = (6 * temp_zone + 177) * M_PI / 180.0;
2761 *Zone = temp_zone;
2762 if (Latitude < 0)
2763 *Hemisphere = 'S';
2764 else
2765 *Hemisphere = 'N';
2766 }
2767 } /* END OF if (!Error_Code) */
2768 return (Error_Code);
2769} /* MAG_GetUTMParameters */
2770
2771int MAG_isNaN(double d) { return d != d; }
2772
2773int MAG_RotateMagneticVector(MAGtype_CoordSpherical CoordSpherical,
2774 MAGtype_CoordGeodetic CoordGeodetic,
2775 MAGtype_MagneticResults MagneticResultsSph,
2776 MAGtype_MagneticResults *MagneticResultsGeo)
2777/* Rotate the Magnetic Vectors to Geodetic Coordinates
2778Manoj Nair, June, 2009 Manoj.C.Nair@Noaa.Gov
2779Equation 16, WMM Technical report
2780
2781INPUT : CoordSpherical : Data structure MAGtype_CoordSpherical with the
2782following elements double lambda; ( longitude) double phig; ( geocentric
2783latitude ) double r; ( distance from the center of the ellipsoid)
2784
2785 CoordGeodetic : Data structure MAGtype_CoordGeodetic with the
2786following elements double lambda; (longitude) double phi; ( geodetic latitude)
2787 double HeightAboveEllipsoid; (height above the ellipsoid
2788(HaE) ) double HeightAboveGeoid;(height above the Geoid )
2789
2790 MagneticResultsSph : Data structure MAGtype_MagneticResults with
2791the following elements double Bx; North double By; East double Bz;
2792Down
2793
2794OUTPUT: MagneticResultsGeo Pointer to the data structure
2795MAGtype_MagneticResults, with the following elements double Bx; North
2796 double By; East
2797 double Bz; Down
2798
2799CALLS : none
2800
2801 */
2802{
2803 double Psi;
2804 /* Difference between the spherical and Geodetic latitudes */
2805 Psi = (M_PI / 180) * (CoordSpherical.phig - CoordGeodetic.phi);
2806
2807 /* Rotate spherical field components to the Geodetic system */
2808 MagneticResultsGeo->Bz =
2809 MagneticResultsSph.Bx * sin(Psi) + MagneticResultsSph.Bz * cos(Psi);
2810 MagneticResultsGeo->Bx =
2811 MagneticResultsSph.Bx * cos(Psi) - MagneticResultsSph.Bz * sin(Psi);
2812 MagneticResultsGeo->By = MagneticResultsSph.By;
2813 return TRUE;
2814} /*MAG_RotateMagneticVector*/
2815
2816void MAG_SphericalToCartesian(MAGtype_CoordSpherical CoordSpherical, double *x,
2817 double *y, double *z) {
2818 double radphi;
2819 double radlambda;
2820
2821 radphi = CoordSpherical.phig * (M_PI / 180);
2822 radlambda = CoordSpherical.lambda * (M_PI / 180);
2823
2824 *x = CoordSpherical.r * cos(radphi) * cos(radlambda);
2825 *y = CoordSpherical.r * cos(radphi) * sin(radlambda);
2826 *z = CoordSpherical.r * sin(radphi);
2827 return;
2828}
2829
2830void MAG_SphericalToGeodetic(MAGtype_Ellipsoid Ellip,
2831 MAGtype_CoordSpherical CoordSpherical,
2832 MAGtype_CoordGeodetic *CoordGeodetic) {
2833 /*This converts spherical coordinates back to geodetic coordinates. It is not
2834 used in the WMM but may be necessary for some applications, such as
2835 geomagnetic coordinates*/
2836 double x, y, z;
2837
2838 MAG_SphericalToCartesian(CoordSpherical, &x, &y, &z);
2839 MAG_CartesianToGeodetic(Ellip, x, y, z, CoordGeodetic);
2840}
2841
2842void MAG_TMfwd4(double Eps, double Epssq, double K0R4, double K0R4oa,
2843 double Acoeff[], double Lam0, double K0, double falseE,
2844 double falseN, int XYonly, double Lambda, double Phi, double *X,
2845 double *Y, double *pscale, double *CoM) {
2846 /* Transverse Mercator forward equations including point-scale and CoM
2847 =--------- =------- =--=--= ---------
2848
2849 Algorithm developed by: C. Rollins August 7, 2006
2850 C software written by: K. Robins
2851
2852
2853 Constants fixed by choice of ellipsoid and choice of projection
2854 parameters
2855 ---------------
2856
2857 Eps Eccentricity (epsilon) of the ellipsoid
2858 Epssq Eccentricity squared
2859 ( R4 Meridional isoperimetric radius )
2860 ( K0 Central scale factor )
2861 K0R4 K0 times R4
2862 K0R4oa K0 times Ratio of R4 over semi-major axis
2863 Acoeff Trig series coefficients, omega as a function of chi
2864 Lam0 Longitude of the central meridian in radians
2865 K0 Central scale factor, for example, 0.9996 for UTM
2866 falseE False easting, for example, 500000 for UTM
2867 falseN False northing
2868
2869 Processing option
2870 ---------- ------
2871
2872 XYonly If one (1), then only X and Y will be properly
2873 computed. Values returned for point-scale
2874 and CoM will merely be the trivial values
2875 for points on the central meridian
2876
2877 Input items that identify the point to be converted
2878 ----- -----
2879
2880 Lambda Longitude (from Greenwich) in radians
2881 Phi Latitude in radians
2882
2883 Output items
2884 ------ -----
2885
2886 X X coordinate (Easting) in meters
2887 Y Y coordinate (Northing) in meters
2888 pscale point-scale (dimensionless)
2889 CoM Convergence-of-meridians in radians
2890 */
2891
2892 double Lam, CLam, SLam, CPhi, SPhi;
2893 double P, part1, part2, denom, CChi, SChi;
2894 double U, V;
2895 double T, Tsq, denom2;
2896 double c2u, s2u, c4u, s4u, c6u, s6u, c8u, s8u;
2897 double c2v, s2v, c4v, s4v, c6v, s6v, c8v, s8v;
2898 double Xstar, Ystar;
2899 double sig1, sig2, comroo;
2900
2901 /*
2902 Ellipsoid to sphere
2903 --------- -- ------
2904
2905 Convert longitude (Greenwhich) to longitude from the central meridian
2906 It is unnecessary to find the (-Pi, Pi] equivalent of the result.
2907 Compute its cosine and sine.
2908 */
2909
2910 Lam = Lambda - Lam0;
2911 CLam = cos(Lam);
2912 SLam = sin(Lam);
2913
2914 /* Latitude */
2915
2916 CPhi = cos(Phi);
2917 SPhi = sin(Phi);
2918
2919 /* Convert geodetic latitude, Phi, to conformal latitude, Chi
2920 Only the cosine and sine of Chi are actually needed. */
2921
2922 P = exp(Eps * ATanH(Eps * SPhi));
2923 part1 = (1 + SPhi) / P;
2924 part2 = (1 - SPhi) * P;
2925 denom = 1 / (part1 + part2);
2926 CChi = 2 * CPhi * denom;
2927 SChi = (part1 - part2) * denom;
2928
2929 /*
2930 Sphere to first plane
2931 ------ -- ----- -----
2932
2933 Apply spherical theory of transverse Mercator to get (u,v) coordinates
2934 Note the order of the arguments in Fortran's version of ArcTan, i.e.
2935 atan2(y, x) = ATan(y/x)
2936 The two argument form of ArcTan is needed here.
2937 */
2938
2939 T = CChi * SLam;
2940 U = ATanH(T);
2941 V = atan2(SChi, CChi * CLam);
2942
2943 /*
2944 Trigonometric multiple angles
2945 ------------- -------- ------
2946
2947 Compute Cosh of even multiples of U
2948 Compute Sinh of even multiples of U
2949 Compute Cos of even multiples of V
2950 Compute Sin of even multiples of V
2951 */
2952
2953 Tsq = T * T;
2954 denom2 = 1 / (1 - Tsq);
2955 c2u = (1 + Tsq) * denom2;
2956 s2u = 2 * T * denom2;
2957 c2v = (-1 + CChi * CChi * (1 + CLam * CLam)) * denom2;
2958 s2v = 2 * CLam * CChi * SChi * denom2;
2959
2960 c4u = 1 + 2 * s2u * s2u;
2961 s4u = 2 * c2u * s2u;
2962 c4v = 1 - 2 * s2v * s2v;
2963 s4v = 2 * c2v * s2v;
2964
2965 c6u = c4u * c2u + s4u * s2u;
2966 s6u = s4u * c2u + c4u * s2u;
2967 c6v = c4v * c2v - s4v * s2v;
2968 s6v = s4v * c2v + c4v * s2v;
2969
2970 c8u = 1 + 2 * s4u * s4u;
2971 s8u = 2 * c4u * s4u;
2972 c8v = 1 - 2 * s4v * s4v;
2973 s8v = 2 * c4v * s4v;
2974
2975 /* First plane to second plane
2976 ----- ----- -- ------ -----
2977
2978 Accumulate terms for X and Y
2979 */
2980
2981 Xstar = Acoeff[3] * s8u * c8v;
2982 Xstar = Xstar + Acoeff[2] * s6u * c6v;
2983 Xstar = Xstar + Acoeff[1] * s4u * c4v;
2984 Xstar = Xstar + Acoeff[0] * s2u * c2v;
2985 Xstar = Xstar + U;
2986
2987 Ystar = Acoeff[3] * c8u * s8v;
2988 Ystar = Ystar + Acoeff[2] * c6u * s6v;
2989 Ystar = Ystar + Acoeff[1] * c4u * s4v;
2990 Ystar = Ystar + Acoeff[0] * c2u * s2v;
2991 Ystar = Ystar + V;
2992
2993 /* Apply isoperimetric radius, scale adjustment, and offsets */
2994
2995 *X = K0R4 * Xstar + falseE;
2996 *Y = K0R4 * Ystar + falseN;
2997
2998 /* Point-scale and CoM
2999 ----- ----- --- --- */
3000
3001 if (XYonly == 1) {
3002 *pscale = K0;
3003 *CoM = 0;
3004 } else {
3005 sig1 = 8 * Acoeff[3] * c8u * c8v;
3006 sig1 = sig1 + 6 * Acoeff[2] * c6u * c6v;
3007 sig1 = sig1 + 4 * Acoeff[1] * c4u * c4v;
3008 sig1 = sig1 + 2 * Acoeff[0] * c2u * c2v;
3009 sig1 = sig1 + 1;
3010
3011 sig2 = 8 * Acoeff[3] * s8u * s8v;
3012 sig2 = sig2 + 6 * Acoeff[2] * s6u * s6v;
3013 sig2 = sig2 + 4 * Acoeff[1] * s4u * s4v;
3014 sig2 = sig2 + 2 * Acoeff[0] * s2u * s2v;
3015
3016 /* Combined square roots */
3017 comroo =
3018 sqrt((1 - Epssq * SPhi * SPhi) * denom2 * (sig1 * sig1 + sig2 * sig2));
3019
3020 *pscale = K0R4oa * 2 * denom * comroo;
3021 *CoM = atan2(SChi * SLam, CLam) + atan2(sig2, sig1);
3022 }
3023} /*MAG_TMfwd4*/
3024
3025#ifndef OPENCPN
3026int MAG_YearToDate(MAGtype_Date *CalendarDate)
3027
3028/* Converts a given Decimal year into a Year, Month and Date
3029it also outputs an error string if there is a problem
3030INPUT CalendarDate Pointer to the data structure with the following elements
3031 double DecimalYear; decimal years
3032OUTPUT CalendarDate Pointer to the data structure with the following
3033elements updated
3034 * int Year
3035 * int Month
3036 * int Day
3037 Error pointer to an error string
3038CALLS : none
3039
3040 */
3041{
3042 int MonthDays[13], CumulativeDays = 0;
3043 int ExtraDay = 0;
3044 int i, DayOfTheYear;
3045
3046 if (CalendarDate->DecimalYear == 0) {
3047 CalendarDate->Year = 0;
3048 CalendarDate->Month = 0;
3049 CalendarDate->Day = 0;
3050 return FALSE;
3051 }
3052
3053 CalendarDate->Year = (int)floor(CalendarDate->DecimalYear);
3054
3055 if ((CalendarDate->Year % 4 == 0 && CalendarDate->Year % 100 != 0) ||
3056 CalendarDate->Year % 400 == 0)
3057 ExtraDay = 1;
3058
3059 DayOfTheYear =
3060 floor((CalendarDate->DecimalYear - (double)CalendarDate->Year) *
3061 (365.0 + (double)ExtraDay) +
3062 0.5) +
3063 1;
3064 /*The above floor is used for rounding, this only works for positive
3065 * integers*/
3066
3067 MonthDays[0] = 0;
3068 MonthDays[1] = 31;
3069 MonthDays[2] = 28 + ExtraDay;
3070 MonthDays[3] = 31;
3071 MonthDays[4] = 30;
3072 MonthDays[5] = 31;
3073 MonthDays[6] = 30;
3074 MonthDays[7] = 31;
3075 MonthDays[8] = 31;
3076 MonthDays[9] = 30;
3077 MonthDays[10] = 31;
3078 MonthDays[11] = 30;
3079 MonthDays[12] = 31;
3080
3081 for (i = 1; i <= 12; i++) {
3082 CumulativeDays = CumulativeDays + MonthDays[i];
3083
3084 if (DayOfTheYear <= CumulativeDays) {
3085 CalendarDate->Month = i;
3086 CalendarDate->Day = MonthDays[i] - (CumulativeDays - DayOfTheYear);
3087 break;
3088 }
3089 }
3090
3091 return TRUE;
3092
3093} /*MAG_YearToDate*/
3094
3095#endif /* OPENCPN */
3096
3097/******************************************************************************
3098 ********************************Spherical Harmonics***************************
3099 * This grouping consists of functions that together take gauss coefficients
3100 * and return a magnetic vector for an input location in spherical coordinates
3101 ******************************************************************************/
3102
3103int MAG_AssociatedLegendreFunction(MAGtype_CoordSpherical CoordSpherical,
3104 int nMax,
3105 MAGtype_LegendreFunction *LegendreFunction)
3106
3107/* Computes all of the Schmidt-semi normalized associated Legendre
3108functions up to degree nMax. If nMax <= 16, function MAG_PcupLow is used.
3109Otherwise MAG_PcupHigh is called.
3110INPUT CoordSpherical A data structure with the following elements
3111 double lambda; ( longitude)
3112 double phig; ( geocentric
3113latitude ) double r; ( distance from the center of the ellipsoid) nMax
3114integer ( Maxumum degree of spherical harmonic secular model)
3115 LegendreFunction Pointer to data structure with the following
3116elements double *Pcup; ( pointer to store Legendre Function ) double *dPcup;
3117( pointer to store Derivative of Lagendre function )
3118
3119OUTPUT LegendreFunction Calculated Legendre variables in the data structure
3120
3121 */
3122{
3123 double sin_phi;
3124 int FLAG = 1;
3125
3126 sin_phi = sin(DEG2RAD(CoordSpherical.phig)); /* sin (geocentric latitude) */
3127
3128 if (nMax <= 16 || (1 - fabs(sin_phi)) <
3129 1.0e-10) /* If nMax is less tha 16 or at the poles */
3130 FLAG = MAG_PcupLow(LegendreFunction->Pcup, LegendreFunction->dPcup, sin_phi,
3131 nMax);
3132 else
3133 FLAG = MAG_PcupHigh(LegendreFunction->Pcup, LegendreFunction->dPcup,
3134 sin_phi, nMax);
3135 if (FLAG == 0) /* Error while computing Legendre variables*/
3136 return FALSE;
3137
3138 return TRUE;
3139} /*MAG_AssociatedLegendreFunction */
3140
3141int MAG_CheckGeographicPole(MAGtype_CoordGeodetic *CoordGeodetic)
3142
3143/* Check if the latitude is equal to -90 or 90. If it is,
3144offset it by 1e-5 to avoid division by zero. This is not currently used in the
3145Geomagnetic main function. This may be used to avoid calling
3146MAG_SummationSpecial. The function updates the input data structure.
3147
3148INPUT CoordGeodetic Pointer to the data structure with the following
3149elements double lambda; (longitude) double phi; ( geodetic latitude) double
3150HeightAboveEllipsoid; (height above the ellipsoid (HaE) ) double
3151HeightAboveGeoid;(height above the Geoid ) OUTPUT CoordGeodetic Pointer to the
3152data structure with the following elements updates double phi; ( geodetic
3153latitude) CALLS : none
3154
3155 */
3156{
3157 CoordGeodetic->phi = CoordGeodetic->phi < (-90.0 + MAG_GEO_POLE_TOLERANCE)
3158 ? (-90.0 + MAG_GEO_POLE_TOLERANCE)
3159 : CoordGeodetic->phi;
3160 CoordGeodetic->phi = CoordGeodetic->phi > (90.0 - MAG_GEO_POLE_TOLERANCE)
3161 ? (90.0 - MAG_GEO_POLE_TOLERANCE)
3162 : CoordGeodetic->phi;
3163 return TRUE;
3164} /*MAG_CheckGeographicPole*/
3165
3166int MAG_ComputeSphericalHarmonicVariables(
3167 MAGtype_Ellipsoid Ellip, MAGtype_CoordSpherical CoordSpherical, int nMax,
3168 MAGtype_SphericalHarmonicVariables *SphVariables)
3169
3170/* Computes Spherical variables
3171 Variables computed are (a/r)^(n+2), cos_m(lamda) and sin_m(lambda) for
3172 spherical harmonic summations. (Equations 10-12 in the WMM Technical Report)
3173 INPUT Ellip data structure with the following elements
3174 double a; semi-major axis of the ellipsoid
3175 double b; semi-minor axis of the ellipsoid
3176 double fla; flattening
3177 double epssq; first eccentricity squared
3178 double eps; first eccentricity
3179 double re; mean radius of ellipsoid
3180 CoordSpherical A data structure with the following
3181 elements double lambda; ( longitude) double phig; ( geocentric latitude )
3182 double r; ( distance from the center of the
3183 ellipsoid)
3184 nMax integer ( Maxumum degree of spherical harmonic
3185 secular model)\
3186
3187 OUTPUT SphVariables Pointer to the data structure with the following
3188 elements double RelativeRadiusPower[MAG_MAX_MODEL_DEGREES+1];
3189 [earth_reference_radius_km sph. radius ]^n double
3190 cos_mlambda[MAG_MAX_MODEL_DEGREES+1]; cp(m) - cosine of (mspherical coord.
3191 longitude) double sin_mlambda[MAG_MAX_MODEL_DEGREES+1]; sp(m) - sine of
3192 (mspherical coord. longitude) CALLS : none
3193 */
3194{
3195 double cos_lambda, sin_lambda;
3196 int m, n;
3197 cos_lambda = cos(DEG2RAD(CoordSpherical.lambda));
3198 sin_lambda = sin(DEG2RAD(CoordSpherical.lambda));
3199 /* for n = 0 ... model_order, compute (Radius of Earth / Spherical radius
3200 r)^(n+2) for n 1..nMax-1 (this is much faster than calling pow MAX_N+1
3201 times). */
3202 SphVariables->RelativeRadiusPower[0] =
3203 (Ellip.re / CoordSpherical.r) * (Ellip.re / CoordSpherical.r);
3204 for (n = 1; n <= nMax; n++) {
3205 SphVariables->RelativeRadiusPower[n] =
3206 SphVariables->RelativeRadiusPower[n - 1] *
3207 (Ellip.re / CoordSpherical.r);
3208 }
3209
3210 /*
3211 Compute cos(m*lambda), sin(m*lambda) for m = 0 ... nMax
3212 cos(a + b) = cos(a)*cos(b) - sin(a)*sin(b)
3213 sin(a + b) = cos(a)*sin(b) + sin(a)*cos(b)
3214 */
3215 SphVariables->cos_mlambda[0] = 1.0;
3216 SphVariables->sin_mlambda[0] = 0.0;
3217
3218 SphVariables->cos_mlambda[1] = cos_lambda;
3219 SphVariables->sin_mlambda[1] = sin_lambda;
3220 for (m = 2; m <= nMax; m++) {
3221 SphVariables->cos_mlambda[m] =
3222 SphVariables->cos_mlambda[m - 1] * cos_lambda -
3223 SphVariables->sin_mlambda[m - 1] * sin_lambda;
3224 SphVariables->sin_mlambda[m] =
3225 SphVariables->cos_mlambda[m - 1] * sin_lambda +
3226 SphVariables->sin_mlambda[m - 1] * cos_lambda;
3227 }
3228 return TRUE;
3229} /*MAG_ComputeSphericalHarmonicVariables*/
3230
3231void MAG_GradY(MAGtype_Ellipsoid Ellip, MAGtype_CoordSpherical CoordSpherical,
3232 MAGtype_CoordGeodetic CoordGeodetic,
3233 MAGtype_MagneticModel *TimedMagneticModel,
3234 MAGtype_GeoMagneticElements GeoMagneticElements,
3235 MAGtype_GeoMagneticElements *GradYElements) {
3236 MAGtype_LegendreFunction *LegendreFunction;
3237 MAGtype_SphericalHarmonicVariables *SphVariables;
3238 int NumTerms;
3239 MAGtype_MagneticResults GradYResultsSph, GradYResultsGeo;
3240
3241 NumTerms =
3242 ((TimedMagneticModel->nMax + 1) * (TimedMagneticModel->nMax + 2) / 2);
3243 LegendreFunction = MAG_AllocateLegendreFunctionMemory(
3244 NumTerms); /* For storing the ALF functions */
3245 SphVariables = MAG_AllocateSphVarMemory(TimedMagneticModel->nMax);
3246 MAG_ComputeSphericalHarmonicVariables(
3247 Ellip, CoordSpherical, TimedMagneticModel->nMax,
3248 SphVariables); /* Compute Spherical Harmonic variables */
3249 MAG_AssociatedLegendreFunction(CoordSpherical, TimedMagneticModel->nMax,
3250 LegendreFunction); /* Compute ALF */
3251 MAG_GradYSummation(
3252 LegendreFunction, TimedMagneticModel, *SphVariables, CoordSpherical,
3253 &GradYResultsSph); /* Accumulate the spherical harmonic coefficients*/
3254 MAG_RotateMagneticVector(
3255 CoordSpherical, CoordGeodetic, GradYResultsSph,
3256 &GradYResultsGeo); /* Map the computed Magnetic fields to Geodetic
3257 coordinates */
3258 MAG_CalculateGradientElements(
3259 GradYResultsGeo, GeoMagneticElements,
3260 GradYElements); /* Calculate the Geomagnetic elements, Equation 18 , WMM
3261 Technical report */
3262
3263 MAG_FreeLegendreMemory(LegendreFunction);
3264 MAG_FreeSphVarMemory(SphVariables);
3265}
3266
3267void MAG_GradYSummation(MAGtype_LegendreFunction *LegendreFunction,
3268 MAGtype_MagneticModel *MagneticModel,
3269 MAGtype_SphericalHarmonicVariables SphVariables,
3270 MAGtype_CoordSpherical CoordSpherical,
3271 MAGtype_MagneticResults *GradY) {
3272 int m, n, index;
3273 double cos_phi;
3274 GradY->Bz = 0.0;
3275 GradY->By = 0.0;
3276 GradY->Bx = 0.0;
3277 for (n = 1; n <= MagneticModel->nMax; n++) {
3278 for (m = 0; m <= n; m++) {
3279 index = (n * (n + 1) / 2 + m);
3280
3281 GradY->Bz -= SphVariables.RelativeRadiusPower[n] *
3282 (-1 * MagneticModel->Main_Field_Coeff_G[index] *
3283 SphVariables.sin_mlambda[m] +
3284 MagneticModel->Main_Field_Coeff_H[index] *
3285 SphVariables.cos_mlambda[m]) *
3286 (double)(n + 1) * (double)(m)*LegendreFunction->Pcup[index] *
3287 (1 / CoordSpherical.r);
3288 GradY->By += SphVariables.RelativeRadiusPower[n] *
3289 (MagneticModel->Main_Field_Coeff_G[index] *
3290 SphVariables.cos_mlambda[m] +
3291 MagneticModel->Main_Field_Coeff_H[index] *
3292 SphVariables.sin_mlambda[m]) *
3293 (double)(m * m) * LegendreFunction->Pcup[index] *
3294 (1 / CoordSpherical.r);
3295 GradY->Bx -= SphVariables.RelativeRadiusPower[n] *
3296 (-1 * MagneticModel->Main_Field_Coeff_G[index] *
3297 SphVariables.sin_mlambda[m] +
3298 MagneticModel->Main_Field_Coeff_H[index] *
3299 SphVariables.cos_mlambda[m]) *
3300 (double)(m)*LegendreFunction->dPcup[index] *
3301 (1 / CoordSpherical.r);
3302 }
3303 }
3304
3305 cos_phi = cos(DEG2RAD(CoordSpherical.phig));
3306 if (fabs(cos_phi) > 1.0e-10) {
3307 GradY->By = GradY->By / (cos_phi * cos_phi);
3308 GradY->Bx = GradY->Bx / (cos_phi);
3309 GradY->Bz = GradY->Bz / (cos_phi);
3310 } else
3311 /* Special calculation for component - By - at Geographic poles.
3312 * If the user wants to avoid using this function, please make sure that
3313 * the latitude is not exactly +/-90. An option is to make use the function
3314 * MAG_CheckGeographicPoles.
3315 */
3316 {
3317 /* MAG_SummationSpecial(MagneticModel, SphVariables, CoordSpherical, GradY);
3318 */
3319 }
3320}
3321
3322int MAG_PcupHigh(double *Pcup, double *dPcup, double x, int nMax)
3323
3324/* This function evaluates all of the Schmidt-semi normalized associated
3325 Legendre functions up to degree nMax. The functions are initially scaled by
3326 10^280 sin^m in order to minimize the effects of underflow at large m
3327 near the poles (see Holmes and Featherstone 2002, J. Geodesy, 76,
3328 279-299). Note that this function performs the same operation as MAG_PcupLow.
3329 However this function also can be used for high degree (large nMax)
3330 models.
3331
3332 Calling Parameters:
3333 INPUT
3334 nMax: Maximum spherical harmonic degree to compute.
3335 x: cos(colatitude) or sin(latitude).
3336
3337 OUTPUT
3338 Pcup: A vector of all associated Legendgre polynomials
3339 evaluated at x up to nMax. The lenght must by greater or equal to
3340 (nMax+1)*(nMax+2)/2. dPcup: Derivative of Pcup(x) with respect to latitude
3341
3342 CALLS : none
3343 Notes:
3344
3345
3346
3347 Adopted from the FORTRAN code written by Mark Wieczorek September 25, 2005.
3348
3349 Manoj Nair, Nov, 2009 Manoj.C.Nair@Noaa.Gov
3350
3351 Change from the previous version
3352 The prevous version computes the derivatives as
3353 dP(n,m)(x)/dx, where x = sin(latitude) (or cos(colatitude) ).
3354 However, the WMM Geomagnetic routines requires dP(n,m)(x)/dlatitude.
3355 Hence the derivatives are multiplied by sin(latitude).
3356 Removed the options for CS phase and normalizations.
3357
3358 Note: In geomagnetism, the derivatives of ALF are usually found with
3359 respect to the colatitudes. Here the derivatives are found with respect
3360 to the latitude. The difference is a sign reversal for the derivative of
3361 the Associated Legendre Functions.
3362
3363 The derivatives can't be computed for latitude = |90| degrees.
3364 */
3365{
3366 double pm2, pm1, pmm, plm, rescalem, z, scalef;
3367 double *f1, *f2, *PreSqr;
3368 int k, kstart, m, n, NumTerms;
3369
3370 NumTerms = ((nMax + 1) * (nMax + 2) / 2);
3371
3372 if (fabs(x) == 1.0) {
3373 printf("Error in PcupHigh: derivative cannot be calculated at poles\n");
3374 return FALSE;
3375 }
3376
3377 f1 = (double *)malloc((NumTerms + 1) * sizeof(double));
3378 if (f1 == NULL) {
3379 MAG_Error(18);
3380 return FALSE;
3381 }
3382
3383 PreSqr = (double *)malloc((NumTerms + 1) * sizeof(double));
3384
3385 if (PreSqr == NULL) {
3386 MAG_Error(18);
3387 return FALSE;
3388 }
3389
3390 f2 = (double *)malloc((NumTerms + 1) * sizeof(double));
3391
3392 if (f2 == NULL) {
3393 MAG_Error(18);
3394 return FALSE;
3395 }
3396
3397 scalef = 1.0e-280;
3398
3399 for (n = 0; n <= 2 * nMax + 1; ++n) {
3400 PreSqr[n] = sqrt((double)(n));
3401 }
3402
3403 k = 2;
3404
3405 for (n = 2; n <= nMax; n++) {
3406 k = k + 1;
3407 f1[k] = (double)(2 * n - 1) / (double)(n);
3408 f2[k] = (double)(n - 1) / (double)(n);
3409 for (m = 1; m <= n - 2; m++) {
3410 k = k + 1;
3411 f1[k] = (double)(2 * n - 1) / PreSqr[n + m] / PreSqr[n - m];
3412 f2[k] =
3413 PreSqr[n - m - 1] * PreSqr[n + m - 1] / PreSqr[n + m] / PreSqr[n - m];
3414 }
3415 k = k + 2;
3416 }
3417
3418 /*z = sin (geocentric latitude) */
3419 z = sqrt((1.0 - x) * (1.0 + x));
3420 pm2 = 1.0;
3421 Pcup[0] = 1.0;
3422 dPcup[0] = 0.0;
3423 if (nMax == 0) return FALSE;
3424 pm1 = x;
3425 Pcup[1] = pm1;
3426 dPcup[1] = z;
3427 k = 1;
3428
3429 for (n = 2; n <= nMax; n++) {
3430 k = k + n;
3431 plm = f1[k] * x * pm1 - f2[k] * pm2;
3432 Pcup[k] = plm;
3433 dPcup[k] = (double)(n) * (pm1 - x * plm) / z;
3434 pm2 = pm1;
3435 pm1 = plm;
3436 }
3437
3438 pmm = PreSqr[2] * scalef;
3439 rescalem = 1.0 / scalef;
3440 kstart = 0;
3441
3442 for (m = 1; m <= nMax - 1; ++m) {
3443 rescalem = rescalem * z;
3444
3445 /* Calculate Pcup(m,m)*/
3446 kstart = kstart + m + 1;
3447 pmm = pmm * PreSqr[2 * m + 1] / PreSqr[2 * m];
3448 Pcup[kstart] = pmm * rescalem / PreSqr[2 * m + 1];
3449 dPcup[kstart] = -((double)(m)*x * Pcup[kstart] / z);
3450 pm2 = pmm / PreSqr[2 * m + 1];
3451 /* Calculate Pcup(m+1,m)*/
3452 k = kstart + m + 1;
3453 pm1 = x * PreSqr[2 * m + 1] * pm2;
3454 Pcup[k] = pm1 * rescalem;
3455 dPcup[k] =
3456 ((pm2 * rescalem) * PreSqr[2 * m + 1] - x * (double)(m + 1) * Pcup[k]) /
3457 z;
3458 /* Calculate Pcup(n,m)*/
3459 for (n = m + 2; n <= nMax; ++n) {
3460 k = k + n;
3461 plm = x * f1[k] * pm1 - f2[k] * pm2;
3462 Pcup[k] = plm * rescalem;
3463 dPcup[k] = (PreSqr[n + m] * PreSqr[n - m] * (pm1 * rescalem) -
3464 (double)(n)*x * Pcup[k]) /
3465 z;
3466 pm2 = pm1;
3467 pm1 = plm;
3468 }
3469 }
3470
3471 /* Calculate Pcup(nMax,nMax)*/
3472 rescalem = rescalem * z;
3473 kstart = kstart + m + 1;
3474 pmm = pmm / PreSqr[2 * nMax];
3475 Pcup[kstart] = pmm * rescalem;
3476 dPcup[kstart] = -(double)(nMax)*x * Pcup[kstart] / z;
3477 free(f1);
3478 free(PreSqr);
3479 free(f2);
3480
3481 return TRUE;
3482} /* MAG_PcupHigh */
3483
3484int MAG_PcupLow(double *Pcup, double *dPcup, double x, int nMax)
3485
3486/* This function evaluates all of the Schmidt-semi normalized associated
3487 Legendre functions up to degree nMax.
3488
3489 Calling Parameters:
3490 INPUT
3491 nMax: Maximum spherical harmonic degree to compute.
3492 x: cos(colatitude) or sin(latitude).
3493
3494 OUTPUT
3495 Pcup: A vector of all associated Legendgre polynomials
3496 evaluated at x up to nMax. dPcup: Derivative of Pcup(x) with respect to
3497 latitude
3498
3499 Notes: Overflow may occur if nMax > 20 , especially for high-latitudes.
3500 Use MAG_PcupHigh for large nMax.
3501
3502 Written by Manoj Nair, June, 2009 . Manoj.C.Nair@Noaa.Gov.
3503
3504 Note: In geomagnetism, the derivatives of ALF are usually found with
3505 respect to the colatitudes. Here the derivatives are found with respect
3506 to the latitude. The difference is a sign reversal for the derivative of
3507 the Associated Legendre Functions.
3508 */
3509{
3510 int n, m, index, index1, index2, NumTerms;
3511 double k, z, *schmidtQuasiNorm;
3512 Pcup[0] = 1.0;
3513 dPcup[0] = 0.0;
3514 /*sin (geocentric latitude) - sin_phi */
3515 z = sqrt((1.0 - x) * (1.0 + x));
3516
3517 NumTerms = ((nMax + 1) * (nMax + 2) / 2);
3518 schmidtQuasiNorm = (double *)malloc((NumTerms + 1) * sizeof(double));
3519
3520 if (schmidtQuasiNorm == NULL) {
3521 MAG_Error(19);
3522 return FALSE;
3523 }
3524
3525 /* First, Compute the Gauss-normalized associated Legendre functions*/
3526 for (n = 1; n <= nMax; n++) {
3527 for (m = 0; m <= n; m++) {
3528 index = (n * (n + 1) / 2 + m);
3529 if (n == m) {
3530 index1 = (n - 1) * n / 2 + m - 1;
3531 Pcup[index] = z * Pcup[index1];
3532 dPcup[index] = z * dPcup[index1] + x * Pcup[index1];
3533 } else if (n == 1 && m == 0) {
3534 index1 = (n - 1) * n / 2 + m;
3535 Pcup[index] = x * Pcup[index1];
3536 dPcup[index] = x * dPcup[index1] - z * Pcup[index1];
3537 } else if (n > 1 && n != m) {
3538 index1 = (n - 2) * (n - 1) / 2 + m;
3539 index2 = (n - 1) * n / 2 + m;
3540 if (m > n - 2) {
3541 Pcup[index] = x * Pcup[index2];
3542 dPcup[index] = x * dPcup[index2] - z * Pcup[index2];
3543 } else {
3544 k = (double)(((n - 1) * (n - 1)) - (m * m)) /
3545 (double)((2 * n - 1) * (2 * n - 3));
3546 Pcup[index] = x * Pcup[index2] - k * Pcup[index1];
3547 dPcup[index] =
3548 x * dPcup[index2] - z * Pcup[index2] - k * dPcup[index1];
3549 }
3550 }
3551 }
3552 }
3553 /* Compute the ration between the the Schmidt quasi-normalized associated
3554 * Legendre functions and the Gauss-normalized version. */
3555
3556 schmidtQuasiNorm[0] = 1.0;
3557 for (n = 1; n <= nMax; n++) {
3558 index = (n * (n + 1) / 2);
3559 index1 = (n - 1) * n / 2;
3560 /* for m = 0 */
3561 schmidtQuasiNorm[index] =
3562 schmidtQuasiNorm[index1] * (double)(2 * n - 1) / (double)n;
3563
3564 for (m = 1; m <= n; m++) {
3565 index = (n * (n + 1) / 2 + m);
3566 index1 = (n * (n + 1) / 2 + m - 1);
3567 schmidtQuasiNorm[index] =
3568 schmidtQuasiNorm[index1] *
3569 sqrt((double)((n - m + 1) * (m == 1 ? 2 : 1)) / (double)(n + m));
3570 }
3571 }
3572
3573 /* Converts the Gauss-normalized associated Legendre
3574 functions to the Schmidt quasi-normalized version using pre-computed
3575 relation stored in the variable schmidtQuasiNorm */
3576
3577 for (n = 1; n <= nMax; n++) {
3578 for (m = 0; m <= n; m++) {
3579 index = (n * (n + 1) / 2 + m);
3580 Pcup[index] = Pcup[index] * schmidtQuasiNorm[index];
3581 dPcup[index] = -dPcup[index] * schmidtQuasiNorm[index];
3582 /* The sign is changed since the new WMM routines use derivative with
3583 respect to latitude insted of co-latitude */
3584 }
3585 }
3586
3587 if (schmidtQuasiNorm) free(schmidtQuasiNorm);
3588 return TRUE;
3589} /*MAG_PcupLow */
3590
3591int MAG_SecVarSummation(MAGtype_LegendreFunction *LegendreFunction,
3592 MAGtype_MagneticModel *MagneticModel,
3593 MAGtype_SphericalHarmonicVariables SphVariables,
3594 MAGtype_CoordSpherical CoordSpherical,
3595 MAGtype_MagneticResults *MagneticResults) {
3596 /*This Function sums the secular variation coefficients to get the secular
3597 variation of the Magnetic vector. INPUT : LegendreFunction MagneticModel
3598 SphVariables
3599 CoordSpherical
3600 OUTPUT : MagneticResults
3601
3602 CALLS : MAG_SecVarSummationSpecial
3603
3604 */
3605 int m, n, index;
3606 double cos_phi;
3607 MagneticModel->SecularVariationUsed = TRUE;
3608 MagneticResults->Bz = 0.0;
3609 MagneticResults->By = 0.0;
3610 MagneticResults->Bx = 0.0;
3611 for (n = 1; n <= MagneticModel->nMaxSecVar; n++) {
3612 for (m = 0; m <= n; m++) {
3613 index = (n * (n + 1) / 2 + m);
3614
3615 /* nMax (n+2) n m m m
3616 Bz = -SUM (a/r) (n+1) SUM [g cos(m p) + h sin(m p)] P
3617 (sin(phi)) n=1 m=0 n n n */
3618 /* Derivative with respect to radius.*/
3619 MagneticResults->Bz -= SphVariables.RelativeRadiusPower[n] *
3620 (MagneticModel->Secular_Var_Coeff_G[index] *
3621 SphVariables.cos_mlambda[m] +
3622 MagneticModel->Secular_Var_Coeff_H[index] *
3623 SphVariables.sin_mlambda[m]) *
3624 (double)(n + 1) * LegendreFunction->Pcup[index];
3625
3626 /* 1 nMax (n+2) n m m m
3627 By = SUM (a/r) (m) SUM [g cos(m p) + h sin(m p)] dP
3628 (sin(phi)) n=1 m=0 n n n */
3629 /* Derivative with respect to longitude, divided by radius. */
3630 MagneticResults->By += SphVariables.RelativeRadiusPower[n] *
3631 (MagneticModel->Secular_Var_Coeff_G[index] *
3632 SphVariables.sin_mlambda[m] -
3633 MagneticModel->Secular_Var_Coeff_H[index] *
3634 SphVariables.cos_mlambda[m]) *
3635 (double)(m)*LegendreFunction->Pcup[index];
3636 /* nMax (n+2) n m m m
3637 Bx = - SUM (a/r) SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
3638 n=1 m=0 n n n */
3639 /* Derivative with respect to latitude, divided by radius. */
3640
3641 MagneticResults->Bx -= SphVariables.RelativeRadiusPower[n] *
3642 (MagneticModel->Secular_Var_Coeff_G[index] *
3643 SphVariables.cos_mlambda[m] +
3644 MagneticModel->Secular_Var_Coeff_H[index] *
3645 SphVariables.sin_mlambda[m]) *
3646 LegendreFunction->dPcup[index];
3647 }
3648 }
3649 cos_phi = cos(DEG2RAD(CoordSpherical.phig));
3650 if (fabs(cos_phi) > 1.0e-10) {
3651 MagneticResults->By = MagneticResults->By / cos_phi;
3652 } else
3653 /* Special calculation for component By at Geographic poles */
3654 {
3655 MAG_SecVarSummationSpecial(MagneticModel, SphVariables, CoordSpherical,
3656 MagneticResults);
3657 }
3658 return TRUE;
3659} /*MAG_SecVarSummation*/
3660
3661int MAG_SecVarSummationSpecial(MAGtype_MagneticModel *MagneticModel,
3662 MAGtype_SphericalHarmonicVariables SphVariables,
3663 MAGtype_CoordSpherical CoordSpherical,
3664 MAGtype_MagneticResults *MagneticResults) {
3665 /*Special calculation for the secular variation summation at the poles.
3666
3667
3668 INPUT: MagneticModel
3669 SphVariables
3670 CoordSpherical
3671 OUTPUT: MagneticResults
3672 CALLS : none
3673
3674
3675 */
3676 int n, index;
3677 double k, sin_phi, *PcupS, schmidtQuasiNorm1, schmidtQuasiNorm2,
3678 schmidtQuasiNorm3;
3679
3680 PcupS = (double *)malloc((MagneticModel->nMaxSecVar + 1) * sizeof(double));
3681
3682 if (PcupS == NULL) {
3683 MAG_Error(15);
3684 return FALSE;
3685 }
3686
3687 PcupS[0] = 1;
3688 schmidtQuasiNorm1 = 1.0;
3689
3690 MagneticResults->By = 0.0;
3691 sin_phi = sin(DEG2RAD(CoordSpherical.phig));
3692
3693 for (n = 1; n <= MagneticModel->nMaxSecVar; n++) {
3694 index = (n * (n + 1) / 2 + 1);
3695 schmidtQuasiNorm2 = schmidtQuasiNorm1 * (double)(2 * n - 1) / (double)n;
3696 schmidtQuasiNorm3 =
3697 schmidtQuasiNorm2 * sqrt((double)(n * 2) / (double)(n + 1));
3698 schmidtQuasiNorm1 = schmidtQuasiNorm2;
3699 if (n == 1) {
3700 PcupS[n] = PcupS[n - 1];
3701 } else {
3702 k = (double)(((n - 1) * (n - 1)) - 1) /
3703 (double)((2 * n - 1) * (2 * n - 3));
3704 PcupS[n] = sin_phi * PcupS[n - 1] - k * PcupS[n - 2];
3705 }
3706
3707 /* 1 nMax (n+2) n m m m
3708 By = SUM (a/r) (m) SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
3709 n=1 m=0 n n n */
3710 /* Derivative with respect to longitude, divided by radius. */
3711
3712 MagneticResults->By += SphVariables.RelativeRadiusPower[n] *
3713 (MagneticModel->Secular_Var_Coeff_G[index] *
3714 SphVariables.sin_mlambda[1] -
3715 MagneticModel->Secular_Var_Coeff_H[index] *
3716 SphVariables.cos_mlambda[1]) *
3717 PcupS[n] * schmidtQuasiNorm3;
3718 }
3719
3720 if (PcupS) free(PcupS);
3721 return TRUE;
3722} /*SecVarSummationSpecial*/
3723
3724int MAG_Summation(MAGtype_LegendreFunction *LegendreFunction,
3725 MAGtype_MagneticModel *MagneticModel,
3726 MAGtype_SphericalHarmonicVariables SphVariables,
3727 MAGtype_CoordSpherical CoordSpherical,
3728 MAGtype_MagneticResults *MagneticResults) {
3729 /* Computes Geomagnetic Field Elements X, Y and Z in Spherical coordinate
3730 system using spherical harmonic summation.
3731
3732
3733 The vector Magnetic field is given by -grad V, where V is Geomagnetic scalar
3734 potential The gradient in spherical coordinates is given by:
3735
3736 dV ^ 1 dV ^ 1 dV ^
3737 grad V = -- r + - -- t + -------- -- p
3738 dr r dt r sin(t) dp
3739
3740
3741 INPUT : LegendreFunction
3742 MagneticModel
3743 SphVariables
3744 CoordSpherical
3745 OUTPUT : MagneticResults
3746
3747 CALLS : MAG_SummationSpecial
3748
3749
3750
3751 Manoj Nair, June, 2009 Manoj.C.Nair@Noaa.Gov
3752 */
3753 int m, n, index;
3754 double cos_phi;
3755 MagneticResults->Bz = 0.0;
3756 MagneticResults->By = 0.0;
3757 MagneticResults->Bx = 0.0;
3758 for (n = 1; n <= MagneticModel->nMax; n++) {
3759 for (m = 0; m <= n; m++) {
3760 index = (n * (n + 1) / 2 + m);
3761
3762 /* nMax (n+2) n m m m
3763 Bz = -SUM (a/r) (n+1) SUM [g cos(m p) + h sin(m p)] P
3764 (sin(phi)) n=1 m=0 n n n */
3765 /* Equation 12 in the WMM Technical report. Derivative with respect to
3766 * radius.*/
3767 MagneticResults->Bz -= SphVariables.RelativeRadiusPower[n] *
3768 (MagneticModel->Main_Field_Coeff_G[index] *
3769 SphVariables.cos_mlambda[m] +
3770 MagneticModel->Main_Field_Coeff_H[index] *
3771 SphVariables.sin_mlambda[m]) *
3772 (double)(n + 1) * LegendreFunction->Pcup[index];
3773
3774 /* 1 nMax (n+2) n m m m
3775 By = SUM (a/r) (m) SUM [g cos(m p) + h sin(m p)] dP
3776 (sin(phi)) n=1 m=0 n n n */
3777 /* Equation 11 in the WMM Technical report. Derivative with respect to
3778 * longitude, divided by radius. */
3779 MagneticResults->By += SphVariables.RelativeRadiusPower[n] *
3780 (MagneticModel->Main_Field_Coeff_G[index] *
3781 SphVariables.sin_mlambda[m] -
3782 MagneticModel->Main_Field_Coeff_H[index] *
3783 SphVariables.cos_mlambda[m]) *
3784 (double)(m)*LegendreFunction->Pcup[index];
3785 /* nMax (n+2) n m m m
3786 Bx = - SUM (a/r) SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
3787 n=1 m=0 n n n */
3788 /* Equation 10 in the WMM Technical report. Derivative with respect to
3789 * latitude, divided by radius. */
3790
3791 MagneticResults->Bx -= SphVariables.RelativeRadiusPower[n] *
3792 (MagneticModel->Main_Field_Coeff_G[index] *
3793 SphVariables.cos_mlambda[m] +
3794 MagneticModel->Main_Field_Coeff_H[index] *
3795 SphVariables.sin_mlambda[m]) *
3796 LegendreFunction->dPcup[index];
3797 }
3798 }
3799
3800 cos_phi = cos(DEG2RAD(CoordSpherical.phig));
3801 if (fabs(cos_phi) > 1.0e-10) {
3802 MagneticResults->By = MagneticResults->By / cos_phi;
3803 } else
3804 /* Special calculation for component - By - at Geographic poles.
3805 * If the user wants to avoid using this function, please make sure that
3806 * the latitude is not exactly +/-90. An option is to make use the function
3807 * MAG_CheckGeographicPoles.
3808 */
3809 {
3810 MAG_SummationSpecial(MagneticModel, SphVariables, CoordSpherical,
3811 MagneticResults);
3812 }
3813 return TRUE;
3814} /*MAG_Summation */
3815
3816int MAG_SummationSpecial(MAGtype_MagneticModel *MagneticModel,
3817 MAGtype_SphericalHarmonicVariables SphVariables,
3818 MAGtype_CoordSpherical CoordSpherical,
3819 MAGtype_MagneticResults *MagneticResults)
3820/* Special calculation for the component By at Geographic poles.
3821Manoj Nair, June, 2009 manoj.c.nair@noaa.gov
3822INPUT: MagneticModel
3823 SphVariables
3824 CoordSpherical
3825OUTPUT: MagneticResults
3826CALLS : none
3827See Section 1.4, "SINGULARITIES AT THE GEOGRAPHIC POLES", WMM Technical report
3828
3829 */
3830{
3831 int n, index;
3832 double k, sin_phi, *PcupS, schmidtQuasiNorm1, schmidtQuasiNorm2,
3833 schmidtQuasiNorm3;
3834
3835 PcupS = (double *)malloc((MagneticModel->nMax + 1) * sizeof(double));
3836 if (PcupS == 0) {
3837 MAG_Error(14);
3838 return FALSE;
3839 }
3840
3841 PcupS[0] = 1;
3842 schmidtQuasiNorm1 = 1.0;
3843
3844 MagneticResults->By = 0.0;
3845 sin_phi = sin(DEG2RAD(CoordSpherical.phig));
3846
3847 for (n = 1; n <= MagneticModel->nMax; n++) {
3848 /*Compute the ration between the Gauss-normalized associated Legendre
3849functions and the Schmidt quasi-normalized version. This is equivalent to
3850sqrt((m==0?1:2)*(n-m)!/(n+m!))*(2n-1)!!/(n-m)! */
3851
3852 index = (n * (n + 1) / 2 + 1);
3853 schmidtQuasiNorm2 = schmidtQuasiNorm1 * (double)(2 * n - 1) / (double)n;
3854 schmidtQuasiNorm3 =
3855 schmidtQuasiNorm2 * sqrt((double)(n * 2) / (double)(n + 1));
3856 schmidtQuasiNorm1 = schmidtQuasiNorm2;
3857 if (n == 1) {
3858 PcupS[n] = PcupS[n - 1];
3859 } else {
3860 k = (double)(((n - 1) * (n - 1)) - 1) /
3861 (double)((2 * n - 1) * (2 * n - 3));
3862 PcupS[n] = sin_phi * PcupS[n - 1] - k * PcupS[n - 2];
3863 }
3864
3865 /* 1 nMax (n+2) n m m m
3866 By = SUM (a/r) (m) SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
3867 n=1 m=0 n n n */
3868 /* Equation 11 in the WMM Technical report. Derivative with respect to
3869 * longitude, divided by radius. */
3870
3871 MagneticResults->By += SphVariables.RelativeRadiusPower[n] *
3872 (MagneticModel->Main_Field_Coeff_G[index] *
3873 SphVariables.sin_mlambda[1] -
3874 MagneticModel->Main_Field_Coeff_H[index] *
3875 SphVariables.cos_mlambda[1]) *
3876 PcupS[n] * schmidtQuasiNorm3;
3877 }
3878
3879 if (PcupS) free(PcupS);
3880 return TRUE;
3881} /*MAG_SummationSpecial */
3882
3883int MAG_TimelyModifyMagneticModel(MAGtype_Date UserDate,
3884 MAGtype_MagneticModel *MagneticModel,
3885 MAGtype_MagneticModel *TimedMagneticModel)
3886
3887/* Time change the Model coefficients from the base year of the model using
3888secular variation coefficients. Store the coefficients of the static model with
3889their values advanced from epoch t0 to epoch t. Copy the SV coefficients. If
3890input "t�" is the same as "t0", then this is merely a copy operation. If the
3891address of "TimedMagneticModel" is the same as the address of "MagneticModel",
3892then this procedure overwrites the given item "MagneticModel".
3893
3894INPUT: UserDate
3895 MagneticModel
3896OUTPUT:TimedMagneticModel
3897CALLS : none
3898 */
3899{
3900 int n, m, index, a, b;
3901 TimedMagneticModel->EditionDate = MagneticModel->EditionDate;
3902 TimedMagneticModel->epoch = MagneticModel->epoch;
3903 TimedMagneticModel->nMax = MagneticModel->nMax;
3904 TimedMagneticModel->nMaxSecVar = MagneticModel->nMaxSecVar;
3905 a = TimedMagneticModel->nMaxSecVar;
3906 b = (a * (a + 1) / 2 + a);
3907 strcpy(TimedMagneticModel->ModelName, MagneticModel->ModelName);
3908 for (n = 1; n <= MagneticModel->nMax; n++) {
3909 for (m = 0; m <= n; m++) {
3910 index = (n * (n + 1) / 2 + m);
3911 if (index <= b) {
3912 TimedMagneticModel->Main_Field_Coeff_H[index] =
3913 MagneticModel->Main_Field_Coeff_H[index] +
3914 (UserDate.DecimalYear - MagneticModel->epoch) *
3915 MagneticModel->Secular_Var_Coeff_H[index];
3916 TimedMagneticModel->Main_Field_Coeff_G[index] =
3917 MagneticModel->Main_Field_Coeff_G[index] +
3918 (UserDate.DecimalYear - MagneticModel->epoch) *
3919 MagneticModel->Secular_Var_Coeff_G[index];
3920 TimedMagneticModel->Secular_Var_Coeff_H[index] =
3921 MagneticModel
3922 ->Secular_Var_Coeff_H[index]; /* We need a copy of the secular
3923 var coef to calculate secular
3924 change */
3925 TimedMagneticModel->Secular_Var_Coeff_G[index] =
3926 MagneticModel->Secular_Var_Coeff_G[index];
3927 } else {
3928 TimedMagneticModel->Main_Field_Coeff_H[index] =
3929 MagneticModel->Main_Field_Coeff_H[index];
3930 TimedMagneticModel->Main_Field_Coeff_G[index] =
3931 MagneticModel->Main_Field_Coeff_G[index];
3932 }
3933 }
3934 }
3935 return TRUE;
3936} /* MAG_TimelyModifyMagneticModel */
3937
3938/*End of Spherical Harmonic Functions*/
3939
3940/******************************************************************************
3941 *************************************Geoid************************************
3942 * This grouping consists of functions that make calculations to adjust
3943 * ellipsoid height to height above the geoid (Height above MSL).
3944 ******************************************************************************
3945 ******************************************************************************/
3946
3947int MAG_ConvertGeoidToEllipsoidHeight(MAGtype_CoordGeodetic *CoordGeodetic,
3948 MAGtype_Geoid *Geoid)
3949
3950/*
3951 * The function Convert_Geoid_To_Ellipsoid_Height converts the specified WGS84
3952 * Geoid height at the specified geodetic coordinates to the equivalent
3953 * ellipsoid height, using the EGM96 gravity model.
3954 *
3955 * CoordGeodetic->phi : Geodetic latitude in degress (input)
3956 * CoordGeodetic->lambda : Geodetic longitude in degrees (input)
3957 * CoordGeodetic->HeightAboveEllipsoid : Ellipsoid height, in
3958 kilometers (output)
3959 * CoordGeodetic->HeightAboveGeoid: Geoid height, in kilometers (input)
3960 *
3961 CALLS : MAG_GetGeoidHeight (
3962
3963 */
3964{
3965 double DeltaHeight;
3966 int Error_Code;
3967 double lat, lon;
3968
3969 if (Geoid->UseGeoid == 1) { /* Geoid correction required */
3970 /* To ensure that latitude is less than 90 call MAG_EquivalentLatLon() */
3971 MAG_EquivalentLatLon(CoordGeodetic->phi, CoordGeodetic->lambda, &lat, &lon);
3972 Error_Code = MAG_GetGeoidHeight(lat, lon, &DeltaHeight, Geoid);
3973 CoordGeodetic->HeightAboveEllipsoid =
3974 CoordGeodetic->HeightAboveGeoid +
3975 DeltaHeight / 1000; /* Input and output should be kilometers,
3976However MAG_GetGeoidHeight returns Geoid height in meters - Hence division by
39771000 */
3978 } else /* Geoid correction not required, copy the MSL height to Ellipsoid
3979 height */
3980 {
3981 CoordGeodetic->HeightAboveEllipsoid = CoordGeodetic->HeightAboveGeoid;
3982 Error_Code = TRUE;
3983 }
3984 return (Error_Code);
3985} /* MAG_ConvertGeoidToEllipsoidHeight*/
3986
3987int MAG_GetGeoidHeight(double Latitude, double Longitude, double *DeltaHeight,
3988 MAGtype_Geoid *Geoid)
3989/*
3990 * The function MAG_GetGeoidHeight returns the height of the
3991 * EGM96 geiod above or below the WGS84 ellipsoid,
3992 * at the specified geodetic coordinates,
3993 * using a grid of height adjustments from the EGM96 gravity model.
3994 *
3995 * Latitude : Geodetic latitude in radians (input)
3996 * Longitude : Geodetic longitude in radians (input)
3997 * DeltaHeight : Height Adjustment, in meters. (output)
3998 * Geoid : MAGtype_Geoid with Geoid grid (input)
3999 CALLS : none
4000 */
4001{
4002 long Index;
4003 double DeltaX, DeltaY;
4004 double ElevationSE, ElevationSW, ElevationNE, ElevationNW;
4005 double OffsetX, OffsetY;
4006 double PostX, PostY;
4007 double UpperY, LowerY;
4008 int Error_Code = 0;
4009
4010 if (!Geoid->Geoid_Initialized) {
4011 MAG_Error(5);
4012 return (FALSE);
4013 }
4014 if ((Latitude < -90) || (Latitude > 90)) { /* Latitude out of range */
4015 Error_Code |= 1;
4016 }
4017 if ((Longitude < -180) || (Longitude > 360)) { /* Longitude out of range */
4018 Error_Code |= 1;
4019 }
4020
4021 if (!Error_Code) { /* no errors */
4022 /* Compute X and Y Offsets into Geoid Height Array: */
4023
4024 if (Longitude < 0.0) {
4025 OffsetX = (Longitude + 360.0) * Geoid->ScaleFactor;
4026 } else {
4027 OffsetX = Longitude * Geoid->ScaleFactor;
4028 }
4029 OffsetY = (90.0 - Latitude) * Geoid->ScaleFactor;
4030
4031 /* Find Four Nearest Geoid Height Cells for specified Latitude, Longitude;
4032 */
4033 /* Assumes that (0,0) of Geoid Height Array is at Northwest corner: */
4034
4035 PostX = floor(OffsetX);
4036 if ((PostX + 1) == Geoid->NumbGeoidCols) PostX--;
4037 PostY = floor(OffsetY);
4038 if ((PostY + 1) == Geoid->NumbGeoidRows) PostY--;
4039
4040 Index = (long)(PostY * Geoid->NumbGeoidCols + PostX);
4041 ElevationNW = (double)Geoid->GeoidHeightBuffer[Index];
4042 ElevationNE = (double)Geoid->GeoidHeightBuffer[Index + 1];
4043
4044 Index = (long)((PostY + 1) * Geoid->NumbGeoidCols + PostX);
4045 ElevationSW = (double)Geoid->GeoidHeightBuffer[Index];
4046 ElevationSE = (double)Geoid->GeoidHeightBuffer[Index + 1];
4047
4048 /* Perform Bi-Linear Interpolation to compute Height above Ellipsoid: */
4049
4050 DeltaX = OffsetX - PostX;
4051 DeltaY = OffsetY - PostY;
4052
4053 UpperY = ElevationNW + DeltaX * (ElevationNE - ElevationNW);
4054 LowerY = ElevationSW + DeltaX * (ElevationSE - ElevationSW);
4055
4056 *DeltaHeight = UpperY + DeltaY * (LowerY - UpperY);
4057 } else {
4058 MAG_Error(17);
4059 return (FALSE);
4060 }
4061 return TRUE;
4062} /*MAG_GetGeoidHeight*/
4063
4064void MAG_EquivalentLatLon(double lat, double lon, double *repairedLat,
4065 double *repairedLon)
4066/*This function takes a latitude and longitude that are ordinarily out of range
4067 and gives in range values that are equivalent on the Earth's surface. This is
4068 required to get correct values for the geoid function.*/
4069{
4070 double colat;
4071 colat = 90 - lat;
4072 *repairedLon = lon;
4073 if (colat < 0) colat = -colat;
4074 while (colat > 360) {
4075 colat -= 360;
4076 }
4077 if (colat > 180) {
4078 colat -= 180;
4079 *repairedLon = *repairedLon + 180;
4080 }
4081 *repairedLat = 90 - colat;
4082 if (*repairedLon > 360) *repairedLon -= 360;
4083 if (*repairedLon < -180) *repairedLon += 360;
4084}
4085
4086/*End of Geoid Functions*/
4087
4088/*New Error Functions*/
4089
4090void MAG_WMMErrorCalc(double H, MAGtype_GeoMagneticElements *Uncertainty) {
4091 double decl_variable, decl_constant;
4092 Uncertainty->F = WMM_UNCERTAINTY_F;
4093 Uncertainty->H = WMM_UNCERTAINTY_H;
4094 Uncertainty->X = WMM_UNCERTAINTY_X;
4095 Uncertainty->Z = WMM_UNCERTAINTY_Z;
4096 Uncertainty->Incl = WMM_UNCERTAINTY_I;
4097 Uncertainty->Y = WMM_UNCERTAINTY_Y;
4098 decl_variable = (WMM_UNCERTAINTY_D_COEF / H);
4099 decl_constant = (WMM_UNCERTAINTY_D_OFFSET);
4100 Uncertainty->Decl =
4101 sqrt(decl_constant * decl_constant + decl_variable * decl_variable);
4102 if (Uncertainty->Decl > 180) {
4103 Uncertainty->Decl = 180;
4104 }
4105}
4106
4107void MAG_PrintUserDataWithUncertainty(
4108 MAGtype_GeoMagneticElements GeomagElements,
4109 MAGtype_GeoMagneticElements Errors, MAGtype_CoordGeodetic SpaceInput,
4110 MAGtype_Date TimeInput, MAGtype_MagneticModel *MagneticModel,
4111 MAGtype_Geoid *Geoid) {
4112 char DeclString[100];
4113 char InclString[100];
4114 MAG_DegreeToDMSstring(GeomagElements.Incl, 2, InclString);
4115 if (GeomagElements.H < 6000 && GeomagElements.H > 2000)
4116 MAG_Warnings(1, GeomagElements.H, MagneticModel);
4117 if (GeomagElements.H < 2000) MAG_Warnings(2, GeomagElements.H, MagneticModel);
4118 if (MagneticModel->SecularVariationUsed == TRUE) {
4119 MAG_DegreeToDMSstring(GeomagElements.Decl, 2, DeclString);
4120 printf("\n Results For \n\n");
4121 if (SpaceInput.phi < 0)
4122 printf("Latitude %.2fS\n", -SpaceInput.phi);
4123 else
4124 printf("Latitude %.2fN\n", SpaceInput.phi);
4125 if (SpaceInput.lambda < 0)
4126 printf("Longitude %.2fW\n", -SpaceInput.lambda);
4127 else
4128 printf("Longitude %.2fE\n", SpaceInput.lambda);
4129 if (Geoid->UseGeoid == 1)
4130 printf("Altitude: %.2f Kilometers above mean sea level\n",
4131 SpaceInput.HeightAboveGeoid);
4132 else
4133 printf("Altitude: %.2f Kilometers above the WGS-84 ellipsoid\n",
4134 SpaceInput.HeightAboveEllipsoid);
4135 printf("Date: %.1f\n", TimeInput.DecimalYear);
4136 printf("\n Main Field\t\t\tSecular Change\n");
4137 printf("F = %9.1f +/- %5.1f nT\t\t Fdot = %5.1f\tnT/yr\n",
4138 GeomagElements.F, Errors.F, GeomagElements.Fdot);
4139 printf("H = %9.1f +/- %5.1f nT\t\t Hdot = %5.1f\tnT/yr\n",
4140 GeomagElements.H, Errors.H, GeomagElements.Hdot);
4141 printf("X = %9.1f +/- %5.1f nT\t\t Xdot = %5.1f\tnT/yr\n",
4142 GeomagElements.X, Errors.X, GeomagElements.Xdot);
4143 printf("Y = %9.1f +/- %5.1f nT\t\t Ydot = %5.1f\tnT/yr\n",
4144 GeomagElements.Y, Errors.Y, GeomagElements.Ydot);
4145 printf("Z = %9.1f +/- %5.1f nT\t\t Zdot = %5.1f\tnT/yr\n",
4146 GeomagElements.Z, Errors.Z, GeomagElements.Zdot);
4147 if (GeomagElements.Decl < 0)
4148 printf("Decl =%20s (WEST) +/-%3.0f Min Ddot = %.1f\tMin/yr\n",
4149 DeclString, 60 * Errors.Decl, 60 * GeomagElements.Decldot);
4150 else
4151 printf("Decl =%20s (EAST) +/-%3.0f Min Ddot = %.1f\tMin/yr\n",
4152 DeclString, 60 * Errors.Decl, 60 * GeomagElements.Decldot);
4153 if (GeomagElements.Incl < 0)
4154 printf("Incl =%20s (UP) +/-%3.0f Min Idot = %.1f\tMin/yr\n",
4155 InclString, 60 * Errors.Incl, 60 * GeomagElements.Incldot);
4156 else
4157 printf("Incl =%20s (DOWN) +/-%3.0f Min Idot = %.1f\tMin/yr\n",
4158 InclString, 60 * Errors.Incl, 60 * GeomagElements.Incldot);
4159 } else {
4160 MAG_DegreeToDMSstring(GeomagElements.Decl, 2, DeclString);
4161 printf("\n Results For \n\n");
4162 if (SpaceInput.phi < 0)
4163 printf("Latitude %.2fS\n", -SpaceInput.phi);
4164 else
4165 printf("Latitude %.2fN\n", SpaceInput.phi);
4166 if (SpaceInput.lambda < 0)
4167 printf("Longitude %.2fW\n", -SpaceInput.lambda);
4168 else
4169 printf("Longitude %.2fE\n", SpaceInput.lambda);
4170 if (Geoid->UseGeoid == 1)
4171 printf("Altitude: %.2f Kilometers above MSL\n",
4172 SpaceInput.HeightAboveGeoid);
4173 else
4174 printf("Altitude: %.2f Kilometers above WGS-84 Ellipsoid\n",
4175 SpaceInput.HeightAboveEllipsoid);
4176 printf("Date: %.1f\n", TimeInput.DecimalYear);
4177 printf("\n Main Field\n");
4178 printf("F = %-9.1f +/-%5.1f nT\n", GeomagElements.F, Errors.F);
4179 printf("H = %-9.1f +/-%5.1f nT\n", GeomagElements.H, Errors.H);
4180 printf("X = %-9.1f +/-%5.1f nT\n", GeomagElements.X, Errors.X);
4181 printf("Y = %-9.1f +/-%5.1f nT\n", GeomagElements.Y, Errors.Y);
4182 printf("Z = %-9.1f +/-%5.1f nT\n", GeomagElements.Z, Errors.Z);
4183 if (GeomagElements.Decl < 0)
4184 printf("Decl =%20s (WEST)+/-%4f\n", DeclString, 60 * Errors.Decl);
4185 else
4186 printf("Decl =%20s (EAST)+/-%4f\n", DeclString, 60 * Errors.Decl);
4187 if (GeomagElements.Incl < 0)
4188 printf("Incl =%20s (UP)+/-%4f\n", InclString, 60 * Errors.Incl);
4189 else
4190 printf("Incl =%20s (DOWN)+/-%4f\n", InclString, 60 * Errors.Incl);
4191 }
4192
4193 if (SpaceInput.phi <= -55 || SpaceInput.phi >= 55)
4194 /* Print Grid Variation */
4195 {
4196 MAG_DegreeToDMSstring(GeomagElements.GV, 2, InclString);
4197 printf("\n\n Grid variation =%20s\n", InclString);
4198 }
4199
4200} /*MAG_PrintUserDataWithUncertainty*/
4201
4202#ifndef OPENCPN
4203void MAG_GetDeg(char *Query_String, double *latitude, double bounds[2]) {
4204 /*Gets a degree value from the user using the standard input*/
4205 char buffer[64], Error_Message[255];
4206 int done, i, j;
4207
4208 printf("%s", Query_String);
4209 while (NULL == fgets(buffer, 64, stdin)) {
4210 printf("%s", Query_String);
4211 }
4212 for (i = 0, done = 0, j = 0; i <= 64 && !done; i++) {
4213 if (buffer[i] == '.') {
4214 j = sscanf(buffer, "%lf", latitude);
4215 if (j == 1)
4216 done = 1;
4217 else
4218 done = -1;
4219 }
4220 if (buffer[i] == ',') {
4221 if (MAG_ValidateDMSstring(buffer, bounds[0], bounds[1], Error_Message)) {
4222 MAG_DMSstringToDegree(buffer, latitude);
4223 done = 1;
4224 } else
4225 done = -1;
4226 }
4227 if (buffer[i] == ' ') /* This detects if there is a ' ' somewhere in the
4228 string, if there is the program tries to interpret the input as Degrees
4229 Minutes Seconds.*/
4230 {
4231 if (MAG_ValidateDMSstring(buffer, bounds[0], bounds[1], Error_Message)) {
4232 MAG_DMSstringToDegree(buffer, latitude);
4233 done = 1;
4234 } else
4235 done = -1;
4236 }
4237 if (buffer[i] == '\0' || done == -1) {
4238 if (MAG_ValidateDMSstring(buffer, bounds[0], bounds[1], Error_Message) &&
4239 done != -1) {
4240 sscanf(buffer, "%lf", latitude);
4241 done = 1;
4242 } else {
4243 printf("%s", Error_Message);
4244 strcpy(buffer, "");
4245 printf(
4246 "\nError encountered, please re-enter as '(-)DDD,MM,SS' or in "
4247 "Decimal Degrees DD.ddd:\n");
4248 while (NULL == fgets(buffer, 40, stdin)) {
4249 printf(
4250 "\nError encountered, please re-enter as '(-)DDD,MM,SS' or in "
4251 "Decimal Degrees DD.ddd:\n");
4252 }
4253 i = -1;
4254 done = 0;
4255 }
4256 }
4257 }
4258}
4259
4260#endif /* OPENCPN */
4261
4262int MAG_GetAltitude(char *Query_String, MAGtype_Geoid *Geoid,
4263 MAGtype_CoordGeodetic *coords, int bounds[2],
4264 int AltitudeSetting) {
4265 int done, j, UpBoundOn;
4266 char tmp;
4267 char buffer[64];
4268 double value;
4269 done = 0;
4270 if (bounds[1] != NO_ALT_MAX) {
4271 UpBoundOn = TRUE;
4272 } else {
4273 UpBoundOn = FALSE;
4274 }
4275 printf("%s", Query_String);
4276
4277 while (!done) {
4278 strcpy(buffer, "");
4279 while (NULL == fgets(buffer, 40, stdin)) {
4280 printf("%s", Query_String);
4281 }
4282 j = 0;
4283 if ((AltitudeSetting != MSLON) &&
4284 (buffer[0] == 'e' || buffer[0] == 'E' ||
4285 AltitudeSetting ==
4286 WGS84ON)) /* User entered height above WGS-84 ellipsoid, copy it to
4287 CoordGeodetic->HeightAboveEllipsoid */
4288 {
4289 if (buffer[0] == 'e' || buffer[0] == 'E') {
4290 j = sscanf(buffer, "%c%lf", &tmp, &coords->HeightAboveEllipsoid);
4291 } else {
4292 j = sscanf(buffer, "%lf", &coords->HeightAboveEllipsoid);
4293 }
4294 if (j == 2) j = 1;
4295 Geoid->UseGeoid = 0;
4296 coords->HeightAboveGeoid = coords->HeightAboveEllipsoid;
4297 value = coords->HeightAboveEllipsoid;
4298 } else /* User entered height above MSL, convert it to the height above
4299 WGS-84 ellipsoid */
4300 {
4301 Geoid->UseGeoid = 1;
4302 j = sscanf(buffer, "%lf", &coords->HeightAboveGeoid);
4303 MAG_ConvertGeoidToEllipsoidHeight(coords, Geoid);
4304 value = coords->HeightAboveGeoid;
4305 }
4306 if (j == 1)
4307 done = 1;
4308 else
4309 printf("\nIllegal Format, please re-enter as '(-)HHH.hhh:'\n");
4310 if ((value < bounds[0] || (value > bounds[1] && UpBoundOn)) && done == 1) {
4311 if (UpBoundOn) {
4312 done = 0;
4313 printf(
4314 "\nWarning: The value you have entered of %f km for the elevation "
4315 "is outside of the required range.\n",
4316 value);
4317 printf(" An elevation between %d km and %d km is needed. \n", bounds[0],
4318 bounds[1]);
4319 if (AltitudeSetting == WGS84ON) {
4320 printf(
4321 "Please enter height above WGS-84 Ellipsoid (in kilometers):\n");
4322 } else if (AltitudeSetting == MSLON) {
4323 printf("Please enter height above mean sea level (in kilometers):\n");
4324 } else {
4325 printf(
4326 "Please enter height in kilometers (prepend E for height above "
4327 "WGS-84 Ellipsoid):");
4328 }
4329 } else {
4330 switch (MAG_Warnings(3, value, NULL)) {
4331 case 0:
4332 return USER_GAVE_UP;
4333 case 1:
4334 done = 0;
4335 printf("Please enter height above sea level (in kilometers):\n");
4336 break;
4337 case 2:
4338 break;
4339 }
4340 }
4341 }
4342 }
4343 return 0;
4344}
MAGtype_CoordGeodetic
Definition GeomagnetismHeader.h:153
MAGtype_CoordSpherical
Definition GeomagnetismHeader.h:161
MAGtype_Date
Definition GeomagnetismHeader.h:167
MAGtype_Ellipsoid
Definition GeomagnetismHeader.h:144
MAGtype_GeoMagneticElements
Definition GeomagnetismHeader.h:192
MAGtype_Geoid
Definition GeomagnetismHeader.h:213
MAGtype_Gradient
Definition GeomagnetismHeader.h:224
MAGtype_LegendreFunction
Definition GeomagnetismHeader.h:174
MAGtype_MagneticModel
Definition GeomagnetismHeader.h:124
MAGtype_MagneticResults
Definition GeomagnetismHeader.h:179
MAGtype_SphericalHarmonicVariables
Definition GeomagnetismHeader.h:185
MAGtype_UTMParameters
Definition GeomagnetismHeader.h:236