zhangyiss/DRTP
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
DRTP/geomag70_linux/geomag70.c

2486 lines
97 KiB
C
Raw Permalink Normal View History

initial upload
2024-09-12 18:48:58 +08:00
/****************************************************************************/
/* */
/* NGDC's Geomagnetic Field Modeling software for the IGRF and WMM */
/* */
/****************************************************************************/
/* */
/* Disclaimer: This program has undergone limited testing. It is */
/* being distributed unoffically. The National Geophysical Data */
/* Center does not guarantee it's correctness. */
/* */
/****************************************************************************/
/* */
/* Version 7.0: */
/* - input file format changed to */
/* -- accept new DGRF2005 coeffs with 0.01 nT precision */
/* -- make sure all values are separated by blanks */
/* -- swapped n and m: first is degree, second is order */
/* - new my_isnan function improves portablility */
/* - corrected feet to km conversion factor */
/* - fixed date conversion errors for yyyy,mm,dd format */
/* - fixed lon/lat conversion errors for deg,min,sec format */
/* - simplified leap year identification */
/* - changed comment: units of ddot and idot are arc-min/yr */
/* - added note that this program computes the secular variation as */
/* the 1-year difference, rather than the instantaneous change, */
/* which can be slightly different */
/* - clarified that height is above ellipsoid, not above mean sea level */
/* although the difference is negligible for magnetics */
/* - changed main(argv,argc) to usual definition main(argc,argv) */
/* - corrected rounding of angles close to 60 minutes */
/* Thanks to all who provided bug reports and suggested fixes */
/* */
/* Stefan Maus Jan-25-2010 */
/* */
/****************************************************************************/
/* */
/* Version 6.1: */
/* Included option to read coordinates from a file and output the */
/* results to a new file, repeating the input and adding columns */
/* for the output */
/* Stefan Maus Jan-31-2008 */
/* */
/****************************************************************************/
/* */
/* Version 6.0: */
/* Bug fixes for the interpolation between models. Also added warnings */
/* for declination at low H and corrected behaviour at geogr. poles. */
/* Placed print-out commands into separate routines to facilitate */
/* fine-tuning of the tables */
/* Stefan Maus Aug-24-2004 */
/* */
/****************************************************************************/
/* */
/* This program calculates the geomagnetic field values from */
/* a spherical harmonic model. Inputs required by the user are: */
/* a spherical harmonic model data file, coordinate preference, */
/* altitude, date/range-step, latitude, and longitude. */
/* */
/* Spherical Harmonic */
/* Model Data File : Name of the data file containing the */
/* spherical harmonic coefficients of */
/* the chosen model. The model and path */
/* must be less than PATH chars. */
/* */
/* Coordinate Preference : Geodetic (WGS84 latitude and altitude */
/* above ellipsoid (WGS84), */
/* or geocentric (spherical, altitude */
/* measured from the center of the Earth). */
/* */
/* Altitude : Altitude above ellipsoid (WGS84). The */
/* program asks for altitude above mean */
/* sea level, because the altitude above */
/* ellipsoid is not known to most users. */
/* The resulting error is very small and */
/* negligible for most practical purposes. */
/* If geocentric coordinate preference is */
/* used, then the altitude must be in the */
/* range of 6370.20 km - 6971.20 km as */
/* measured from the center of the earth. */
/* Enter altitude in kilometers, meters, */
/* or feet */
/* */
/* Date : Date, in decimal years, for which to */
/* calculate the values of the magnetic */
/* field. The date must be within the */
/* limits of the model chosen. */
/* */
/* Latitude : Entered in decimal degrees in the */
/* form xxx.xxx. Positive for northern */
/* hemisphere, negative for the southern */
/* hemisphere. */
/* */
/* Longitude : Entered in decimal degrees in the */
/* form xxx.xxx. Positive for eastern */
/* hemisphere, negative for the western */
/* hemisphere. */
/* */
/****************************************************************************/
/* */
/* Subroutines called : degrees_to_decimal,julday,getshc,interpsh, */
/* extrapsh,shval3,dihf,safegets */
/* */
/****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
int my_isnan(double d)
{
return (d != d); /* IEEE: only NaN is not equal to itself */
}
#define NaN log(-1.0)
#define FT2KM (1.0/0.0003048)
#define PI 3.141592654
#define RAD2DEG (180.0/PI)
#ifndef SEEK_SET
#define SEEK_SET 0
#define SEEK_CUR 1
#define SEEK_END 2
#endif
#define IEXT 0
#define FALSE 0
#define TRUE 1 /* constants */
#define RECL 81
#define MAXINBUFF RECL+14
/** Max size of in buffer **/
#define MAXREAD MAXINBUFF-2
/** Max to read 2 less than total size (just to be safe) **/
#define MAXMOD 30
/** Max number of models in a file **/
#define PATH MAXREAD
/** Max path and filename length **/
#define EXT_COEFF1 (double)0
#define EXT_COEFF2 (double)0
#define EXT_COEFF3 (double)0
#define MAXDEG 13
#define MAXCOEFF (MAXDEG*(MAXDEG+2)+1) /* index starts with 1!, (from old Fortran?) */
double gh1[MAXCOEFF];
double gh2[MAXCOEFF];
double gha[MAXCOEFF]; /* Geomag global variables */
double ghb[MAXCOEFF];
double d=0,f=0,h=0,i=0;
double dtemp,ftemp,htemp,itemp;
double x=0,y=0,z=0;
double xtemp,ytemp,ztemp;
FILE *stream = NULL; /* Pointer to specified model data file */
/****************************************************************************/
/* */
/* Program Geomag */
/* */
/****************************************************************************/
/* */
/* This program, originally written in FORTRAN, was developed using */
/* subroutines written by */
/* A. Zunde */
/* USGS, MS 964, Box 25046 Federal Center, Denver, Co. 80225 */
/* and */
/* S.R.C. Malin & D.R. Barraclough */
/* Institute of Geological Sciences, United Kingdom. */
/* */
/* Translated */
/* into C by : Craig H. Shaffer */
/* 29 July, 1988 */
/* */
/* Rewritten by : David Owens */
/* For Susan McLean */
/* */
/* Maintained by: Adam Woods */
/* Contact : geomag.models@noaa.gov */
/* National Geophysical Data Center */
/* World Data Center-A for Solid Earth Geophysics */
/* NOAA, E/GC1, 325 Broadway, */
/* Boulder, CO 80303 */
/* */
/* */
/****************************************************************************/
/* */
/* Some variables used in this program */
/* */
/* Name Type Usage */
/* ------------------------------------------------------------------------ */
/* */
/* a2,b2 Scalar Double Squares of semi-major and semi-minor */
/* axes of the reference spheroid used */
/* for transforming between geodetic or */
/* geocentric coordinates. */
/* */
/* minalt Double array of MAXMOD Minimum height of model. */
/* */
/* altmin Double Minimum height of selected model. */
/* */
/* altmax Double array of MAXMOD Maximum height of model. */
/* */
/* maxalt Double Maximum height of selected model. */
/* */
/* d Scalar Double Declination of the field from the */
/* geographic north (deg). */
/* */
/* sdate Scalar Double start date inputted */
/* */
/* ddot Scalar Double annual rate of change of decl. */
/* (arc-min/yr) */
/* */
/* alt Scalar Double altitude above WGS84 Ellipsoid */
/* */
/* epoch Double array of MAXMOD epoch of model. */
/* */
/* ext Scalar Double Three 1st-degree external coeff. */
/* */
/* latitude Scalar Double Latitude. */
/* */
/* longitude Scalar Double Longitude. */
/* */
/* gh1 Double array Schmidt quasi-normal internal */
/* spherical harmonic coeff. */
/* */
/* gh2 Double array Schmidt quasi-normal internal */
/* spherical harmonic coeff. */
/* */
/* gha Double array Coefficients of resulting model. */
/* */
/* ghb Double array Coefficients of rate of change model.*/
/* */
/* i Scalar Double Inclination (deg). */
/* */
/* idot Scalar Double Rate of change of i (arc-min/yr). */
/* */
/* igdgc Integer Flag for geodetic or geocentric */
/* coordinate choice. */
/* */
/* inbuff Char a of MAXINBUF Input buffer. */
/* */
/* irec_pos Integer array of MAXMOD Record counter for header */
/* */
/* stream Integer File handles for an opened file. */
/* */
/* fileline Integer Current line in file (for errors) */
/* */
/* max1 Integer array of MAXMOD Main field coefficient. */
/* */
/* max2 Integer array of MAXMOD Secular variation coefficient. */
/* */
/* max3 Integer array of MAXMOD Acceleration coefficient. */
/* */
/* mdfile Character array of PATH Model file name. */
/* */
/* minyr Double Min year of all models */
/* */
/* maxyr Double Max year of all models */
/* */
/* yrmax Double array of MAXMOD Max year of model. */
/* */
/* yrmin Double array of MAXMOD Min year of model. */
/* */
/****************************************************************************/
int main(int argc, char**argv)
{
#ifdef MAC
ccommand(argc, argv);
#endif
/* Variable declaration */
/* Control variables */
int again = 1;
int decyears = 3;
int units = 4;
int decdeg = 3;
int range = -1;
int counter = 0;
int warn_H, warn_H_strong, warn_P;
int modelI; /* Which model (Index) */
int nmodel; /* Number of models in file */
int max1[MAXMOD];
int max2[MAXMOD];
int max3[MAXMOD];
int nmax;
int igdgc=3;
int isyear=-1;
int ismonth=-1;
int isday=-1;
int ieyear=-1;
int iemonth=-1;
int ieday=-1;
int ilat_deg=200;
int ilat_min=200;
int ilat_sec=200;
int ilon_deg=200;
int ilon_min=200;
int ilon_sec=200;
int fileline;
long irec_pos[MAXMOD];
int coords_from_file = 0;
int arg_err = 0;
int need_to_read_model = 1;
char mdfile[PATH];
char inbuff[MAXINBUFF];
char model[MAXMOD][9];
char *begin;
char *rest;
char args[7][MAXREAD];
int iarg;
char coord_fname[PATH];
char out_fname[PATH];
FILE *coordfile,*outfile;
int iline=0;
int read_flag;
double epoch[MAXMOD];
double yrmin[MAXMOD];
double yrmax[MAXMOD];
double minyr;
double maxyr;
double altmin[MAXMOD];
double altmax[MAXMOD];
double minalt;
double maxalt;
double alt=-999999;
double sdate=-1;
double step=-1;
double syr;
double edate=-1;
double latitude=200;
double longitude=200;
double ddot;
double fdot;
double hdot;
double idot;
double xdot;
double ydot;
double zdot;
double warn_H_val, warn_H_strong_val;
/* Subroutines used */
void print_dashed_line();
void print_long_dashed_line(void);
void print_header();
void print_result(double date, double d, double i, double h, double x, double y, double z, double f);
void print_header_sv();
void print_result_sv(double date, double ddot, double idot, double hdot, double xdot, double ydot, double zdot, double fdot);
void print_result_file(FILE *outf, double d, double i, double h, double x, double y, double z, double f,
double ddot, double idot, double hdot, double xdot, double ydot, double zdot, double fdot);
double degrees_to_decimal();
double julday();
int interpsh();
int extrapsh();
int shval3();
int dihf();
int safegets(char *buffer,int n);
int getshc();
/* Initializations. */
inbuff[MAXREAD+1]='\0'; /* Just to protect mem. */
inbuff[MAXINBUFF-1]='\0'; /* Just to protect mem. */
for (iarg=0; iarg<argc; iarg++)
if (argv[iarg] != NULL)
strncpy(args[iarg],argv[iarg],MAXREAD);
/* printing out version number and header */
printf("\n\n Geomag v7.0 - Jan 25, 2010 ");
if ((argc==2)&&((*(args[1])=='h')||(*(args[1])=='?')||(args[1][1]=='?')))
{
printf("\n\nUSAGE:\n");
printf("interactive: geomag\n");
printf("command line: geomag model_file date coord alt lat lon\n");
printf("coordinate file: geomag model_file f input_file output_file\n");
printf("or for help: geomag h (or ? or -? or /?) \n");
printf("\n");
printf("Date and location Formats: \n");
printf(" Date: xxxx.xxx for decimal (1947.32)\n");
printf(" YYYY,MM,DD for year, month, day (1947,10,13)\n");
printf(" or start_date-end_date (1943.21-1966.11)\n");
printf(" or start_date-end_date-step (1943.21-1966.11-1.2)\n");
printf(" Coord: D - Geodetic (WGS84 latitude and altitude above mean sea level)\n");
printf(" C - Geocentric (spherical, altitude relative to Earth's center)\n");
printf(" Altitude: Kxxxxxx.xxx for kilometers (K1000.13)\n");
printf(" Mxxxxxx.xxx for meters (m1389.24)\n");
printf(" Fxxxxxx.xxx for feet (F192133.73)\n");
printf(" Lat/Lon: xxx.xxx for decimal (-76.53)\n");
printf(" ddd,mm,ss for degrees, minutes, seconds (-60,23,22)\n");
printf(" (Lat and Lon must be specified in the same format,\n");
printf(" for ddd,mm,ss format, two commas each are required \n");
printf(" and decimals of arc-seconds are ignored) \n");
printf("\nRange (only for interactive and command line options): \n");
printf(" Date and altitude must fit model.\n");
printf(" Lat: -90 to 90 (Use - to denote Southern latitude.)\n");
printf(" Lon: -180 to 180 (Use - to denote Westen longitude.)\n");
printf(" Minutes and Seconds: -59 to 59.\n");
printf("\nNote for command line option: \n");
printf(" All arguments are optional but must preserve order.\n");
printf(" You can only omit arguments AT THE END, not in between.\n");
printf(" Invalid arguments are assumed to be 0.\n");
printf("\nNote that this program computes secular variation as the \n");
printf(" change over the coming year, rather than instantaneous,\n");
printf(" which can be different for declination near magnetic poles.\n");
exit(2);
} /* help */
if ((argc==5)&&(*(args[2])=='f'))
{
printf("\n\n 'f' switch: converting file with multiple locations.\n");
printf(" The first five output columns repeat the input coordinates.\n");
printf(" Then follows D, I, H, X, Y, Z, and F.\n");
printf(" Finally the SV: dD, dI, dH, dX, dY, dZ, and dF\n");
printf(" The units are the same as when the program is\n");
printf(" run in command line or interactive mode.\n\n");
coords_from_file = 1;
strncpy(coord_fname,args[3],MAXREAD);
coordfile=fopen(coord_fname, "rt");
strncpy(out_fname,args[4],MAXREAD);
outfile=fopen(out_fname, "w");
fprintf(outfile,"Date Coord-System Altitude Latitude Longitude D_deg D_min I_deg I_min H_nT X_nT Y_nT Z_nT F_nT dD_min dI_min dH_nT dX_nT dY_nT dZ_nT dF_nT\n");
} /* file option */
if (argc>=3 && argc !=5 && *(args[2])=='f')
{
printf("\n\nERROR in 'f' switch option: wrong number of arguments\n");
}
if ((argc==3)||(argc==4)||(argc==6))
{
printf("\n\nUSAGE:\n");
printf("interactive: geomag\n");
printf("command line: geomag model_file date coord alt lat lon\n");
printf("coordinate file: geomag model_file f input_file output_file\n");
printf("or for help: geomag h \n\n");
} /* missing arguments (not fatal) */
while (again == 1)
{
if (coords_from_file)
{
argc = 7;
read_flag = fscanf(coordfile,"%s%s%s%s%s%*[^\n]",args[2],args[3],args[4],args[5],args[6]);
if (read_flag == EOF) goto reached_EOF;
fprintf(outfile,"%s %s %s %s %s ",args[2],args[3],args[4],args[5],args[6]);fflush(outfile);
iline++;
} /* coords_from_file */
/* Switch on how many arguments are supplied. */
/* Note that there are no 'breaks' after the cases, so these are entry points */
switch(argc)
{
case 7 : strncpy(inbuff, args[6], MAXREAD);
if ((rest=strchr(inbuff, ','))) /* If it contains a comma */
{
decdeg=2; /* Then not decimal degrees */
begin=inbuff;
rest[0]='\0'; /* Chop into sub string */
rest++; /* Move to next substring */
ilon_deg=atoi(begin);
begin=rest;
if ((rest=strchr(begin, ',')))
{
rest[0]='\0';
rest++;
ilon_min=atoi(begin);
ilon_sec=atoi(rest);
}
else
{
ilon_min=0;
ilon_sec=0;
}
}
else
{
decdeg=1; /* Else it's decimal */
longitude=atof(args[6]);
}
case 6 : strncpy(inbuff, args[5], MAXREAD);
if ((rest=strchr(inbuff, ',')))
{
decdeg=2;
begin=inbuff;
rest[0]='\0';
rest++;
ilat_deg=atoi(begin);
begin=rest;
if ((rest=strchr(begin, ',')))
{
rest[0]='\0';
rest++;
ilat_min=atoi(begin);
ilat_sec=atoi(rest);
}
else
{
ilat_min=0;
ilat_sec=0;
}
}
else
{
decdeg=1;
latitude=atof(args[5]);
}
case 5 : strncpy(inbuff, args[4], MAXREAD);
inbuff[0]=toupper(inbuff[0]);
if (inbuff[0]=='K') units=1;
else if (inbuff[0]=='M') units=2;
else if (inbuff[0]=='F') units=3;
if (strlen(inbuff)>1)
{
inbuff[0]='\0';
begin=inbuff+1;
alt=atof(begin);
}
case 4 : strncpy(inbuff, args[3], MAXREAD);
inbuff[0]=toupper(inbuff[0]);
if (inbuff[0]=='D') igdgc=1;
else if (inbuff[0]=='C') igdgc=2;
case 3 : strncpy(inbuff, args[2], MAXREAD);
if ((rest=strchr(inbuff, '-'))) /* If it contains a dash */
{
range = 2; /* They want a range */
rest[0]='\0'; /* Sep dates */
rest++;
begin=rest;
if ((rest=strchr(begin, '-'))) /* If it contains 2 dashs */
{
rest[0]='\0'; /* Sep step */
rest++;
step=atof(rest); /* Get step size */
}
if ((rest=strchr(begin, ','))) /* If it contains a comma */
{
decyears=2; /* It's not decimal years */
rest[0]='\0';
rest++;
ieyear=atoi(begin);
begin=rest;
if ((rest=strchr(begin, ',')))
{
rest[0]='\0';
rest++;
iemonth=atoi(begin);
ieday=atoi(rest);
}
else
{
iemonth=0;
ieday=0;
}
if ((rest=strchr(inbuff, ',')))
{
begin=inbuff;
rest[0]='\0';
rest++;
isyear=atoi(begin);
begin=rest;
if ((rest=strchr(begin, ',')))
{
rest[0]='\0';
rest++;
ismonth=atoi(begin);
isday=atoi(rest);
}
else
{
ismonth=0;
isday=0;
}
}
else
{
sdate=atof(inbuff);
}
}
else
{
decyears=1; /* Else it's decimal years */
sdate=atof(inbuff);
edate=atof(begin);
}
}
else
{
range = 1;
if ((rest=strchr(inbuff, ','))) /* If it contains a comma */
{
decyears=2; /* It's not decimal years */
begin=inbuff;
rest[0]='\0';
rest++;
isyear=atoi(begin);
begin=rest;
if ((rest=strchr(begin, ',')))
{
rest[0]='\0';
rest++;
ismonth=atoi(begin);
isday=atoi(rest);
}
else
{
ismonth=0;
isday=0;
}
sdate = julday(ismonth,isday,isyear);
}
else
{
decyears=1; /* Else it's decimal years */
sdate=atof(args[2]);
}
}
if (sdate==0)
{ /* If date not valid */
decyears=-1;
range=-1;
}
case 2 :
if (need_to_read_model)
{
strncpy(mdfile,args[1],MAXREAD);
stream=fopen(mdfile, "rt");
}
break;
}
if (range == 2 && coords_from_file)
{
printf("Error in line %1d, date = %s: date ranges not allowed for file option\n\n",iline,args[2]);
exit(2);
}
/* Obtain the desired model file and read the data */
warn_H = 0;
warn_H_val = 99999.0;
warn_H_strong = 0;
warn_H_strong_val = 99999.0;
warn_P = 0;
if (need_to_read_model)
{
while (stream == NULL)
{
printf("\n\n");
printf("What is the name of the model data file to be opened? ==> ");
safegets(inbuff,MAXREAD);
strcpy(mdfile, inbuff);
if (!(stream = fopen(mdfile, "rt")))
printf("\nError opening file %s.", mdfile);
}
rewind(stream);
fileline = 0; /* First line will be 1 */
modelI = -1; /* First model will be 0 */
while (fgets(inbuff,MAXREAD,stream)) /* While not end of file
* read to end of line or buffer */
{
fileline++; /* On new line */
if (strlen(inbuff) != RECL) /* IF incorrect record size */
{
printf("Corrupt record in file %s on line %d.\n", mdfile, fileline);
fclose(stream);
exit(5);
}
/* old statement Dec 1999 */
/* if (!strncmp(inbuff," ",4)){ /* If 1st 4 chars are spaces */
/* New statement Dec 1999 changed by wmd required by year 2000 models */
if (!strncmp(inbuff," ",3)) /* If 1st 3 chars are spaces */
{
modelI++; /* New model */
if (modelI > MAXMOD) /* If too many headers */
{
printf("Too many models in file %s on line %d.", mdfile, fileline);
fclose(stream);
exit(6);
}
irec_pos[modelI]=ftell(stream);
/* Get fields from buffer into individual vars. */
sscanf(inbuff, "%s%lg%d%d%d%lg%lg%lg%lg", model[modelI], &epoch[modelI],
&max1[modelI], &max2[modelI], &max3[modelI], &yrmin[modelI],
&yrmax[modelI], &altmin[modelI], &altmax[modelI]);
/* Compute date range for all models */
if (modelI == 0) /*If first model */
{
minyr=yrmin[0];
maxyr=yrmax[0];
}
else
{
if (yrmin[modelI]<minyr)
{
minyr=yrmin[modelI];
}
if (yrmax[modelI]>maxyr){
maxyr=yrmax[modelI];
}
} /* if modelI != 0 */
} /* If 1st 3 chars are spaces */
} /* While not end of model file */
nmodel = modelI + 1;
fclose(stream);
/* if date specified in command line then warn if past end of validity */
if ((sdate>maxyr)&&(sdate<maxyr+1))
{
printf("\nWarning: The date %4.2f is out of range,\n", sdate);
printf(" but still within one year of model expiration date.\n");
printf(" An updated model file is available before 1.1.%4.0f\n",maxyr);
}
} /* if need_to_read_model */
/* Take in field data */
/* Get date */
if (coords_from_file && !arg_err && (decyears != 1 && decyears != 2))
{printf("\nError: unrecognized date %s in coordinate file line %1d\n\n",args[2],iline); arg_err = 1;}
while ((decyears!=1)&&(decyears!=2))
{
printf("\nHow would you like to enter the date?\n");
printf(" 1) In decimal years.\n");
printf(" 2) In year, month, and day.\n");
printf("\n ==> ");
safegets(inbuff, MAXREAD);
decyears=atoi(inbuff);
}
if (coords_from_file && !arg_err && range != 1)
{printf("\nError: unrecognized date %s in coordinate file line %1d\n\n",args[2],iline); arg_err = 1;}
while ((range!=1)&&(range!=2))
{
printf("\nWould you like output for a single date or for a range of dates?\n");
printf(" 1) A single date.\n");
printf(" 2) A range of dates.\n");
printf("\n ==> ");
safegets(inbuff, MAXREAD);
range=atoi(inbuff);
}
if (range == 1)
{
if (coords_from_file && !arg_err && (sdate < minyr || sdate > maxyr+1))
{printf("\nError: unrecognized date %s in coordinate file line %1d\n\n",args[2],iline); arg_err = 1;}
while ((sdate<minyr)||(sdate>maxyr+1))
{
if (decyears==1)
{
printf("\nEnter the decimal date (%4.2f to %4.0f): ",minyr, maxyr);
safegets(inbuff, MAXREAD);
sdate=atof(inbuff);
}
else
{
while ((isyear>(int)maxyr+1)||(isyear<(int)minyr))
{
printf("\nEnter the date (%4.2f to %4.2f)\n ", minyr, maxyr);
printf("\n Year (%d to %d): ",(int)minyr,(int)maxyr);
safegets(inbuff, MAXREAD);
isyear=atoi(inbuff);
}
while ((ismonth>12)||(ismonth<1))
{
printf("\n Month (1-12): ");
safegets(inbuff, MAXREAD);
ismonth=atoi(inbuff);
}
while ((isday>31)||(isday<1))
{
printf("\n Day (1-31): ");
safegets(inbuff, MAXREAD);
isday=atoi(inbuff);
}
sdate = julday(ismonth,isday,isyear);
}
if ((sdate<minyr)||(sdate>=maxyr+1))
{
ismonth=isday=isyear=0;
printf("\nError: The date %4.2f is out of range.\n", sdate);
}
if ((sdate>maxyr)&&(sdate<maxyr+1))
{
printf("\nWarning: The date %4.2f is out of range,\n", sdate);
printf(" but still within one year of model expiration date.\n");
printf(" An updated model file is available before 1.1.%4.0f\n",maxyr);
}
} /* if single date */
} /* (range == 1) */
else
{
while ((sdate<minyr)||(sdate>maxyr))
{
if (decyears==1)
{
printf("\nEnter the decimal start date (%4.2f to %4.0f): ",minyr, maxyr);
safegets(inbuff, MAXREAD);
sdate=atof(inbuff);
}
else
{
while ((isyear>(int)maxyr)||(isyear<(int)minyr))
{
ismonth=isday=isyear=0;
printf("\nEnter the start date (%4.2f to %4.2f)\n ", minyr, maxyr);
printf("\n Year (%d to %d): ",(int)minyr,(int)(maxyr));
safegets(inbuff, MAXREAD);
isyear=atoi(inbuff);
}
while ((ismonth>12)||(ismonth<1))
{
printf("\n Month (1-12): ");
safegets(inbuff, MAXREAD);
ismonth=atoi(inbuff);
}
while ((isday>31)||(isday<1))
{
printf("\n Day (1-31): ");
safegets(inbuff, MAXREAD);
isday=atoi(inbuff);
}
sdate = julday(ismonth,isday,isyear);
if ((sdate<minyr)||(sdate>maxyr))
{
printf("\nThe start date %4.2f is out of range.\n", sdate);
}
}
} /* WHILE ((sdate<minyr)||(sdate>maxyr)) */
while ((edate<=sdate)||(edate>maxyr+1))
{
if (decyears==1)
{
printf("\nEnter the decimal end date (%4.2f to %4.0f): ",sdate, maxyr);
safegets(inbuff, MAXREAD);
edate=atof(inbuff);
}
else
{
while ((ieyear>(int)maxyr)||(ieyear<(int)sdate))
{
iemonth=ieday=ieyear=0;
printf("\nEnter the end date (%4.2f to %4.0f)\n ", sdate, maxyr);
printf("\n Year (%d to %d): ",(int)sdate,(int)(maxyr));
safegets(inbuff, MAXREAD);
ieyear=atoi(inbuff);
}
while ((iemonth>12)||(iemonth<1))
{
printf("\n Month (1-12): ");
safegets(inbuff, MAXREAD);
iemonth=atoi(inbuff);
}
while ((ieday>31)||(ieday<1))
{
printf("\n Day (1-31): ");
safegets(inbuff, MAXREAD);
ieday=atoi(inbuff);
}
edate = julday(iemonth,ieday,ieyear);
}
if ((edate<sdate)||(edate>maxyr+1))
{
printf("\nThe date %4.2f is out of range.\n", edate);
}
if (edate>maxyr && edate<=maxyr+1)
{
printf("\nWarning: The end date %4.2f is out of range,\n", edate);
printf(" but still within one year of model expiration date.\n");
printf(" An updated model file is available before 1.1.%4.0f\n",maxyr);
}
} /* while ((edate<=sdate)||(edate>maxyr+1)) */
while ((step<=0)||(step>(edate-sdate)))
{
printf("\nEnter the step size in years. (0 to %4.2f): ",edate-sdate);
safegets(inbuff, MAXREAD);
step=atof(inbuff);
}
} /* if (range == 2) */
/* Pick model */
for (modelI=0; modelI<nmodel; modelI++)
if (sdate<yrmax[modelI]) break;
if (modelI == nmodel) modelI--; /* if beyond end of last model use last model */
/* Get altitude min and max for selected model. */
minalt=altmin[modelI];
maxalt=altmax[modelI];
/* Get Coordinate prefs */
if (coords_from_file && !arg_err && (igdgc != 1 && igdgc != 2))
{printf("\nError: unrecognized coordinate system %s in coordinate file line %1d\n\n",args[3],iline); arg_err = 1;}
while ((igdgc!=1)&&(igdgc!=2))
{
printf("\n\nEnter Coordinate Preferences:");
printf("\n 1) Geodetic (WGS84 latitude and altitude above mean sea level)");
printf("\n 2) Geocentric (spherical, altitude relative to Earth's center)\n");
printf("\n ==> ");
safegets(inbuff, MAXREAD);
igdgc=atoi(inbuff);
}
/* If needed modify ranges to reflect coords. */
if (igdgc==2)
{
minalt+=6371.2; /* Add radius to ranges. */
maxalt+=6371.2;
}
/* Get unit prefs */
if (igdgc==1)
{
if (coords_from_file && !arg_err && (units > 3 || units < 1))
{printf("\nError: unrecognized altitude units %s in coordinate file line %1d\n\n",args[4],iline); arg_err = 1;}
while ((units>3)||(units<1))
{
printf("\n\nEnter Unit Preferences:");
printf("\n 1) Kilometers");
printf("\n 2) Meters");
printf("\n 3) Feet\n");
printf("\n ==> ");
safegets(inbuff, MAXREAD);
units=atoi(inbuff);
}
}
else units = 1; /* geocentric always in km */
/* Do unit conversions if neccessary */
if (units==2)
{
minalt*=1000.0;
maxalt*=1000.0;
}
else if (units==3)
{
minalt*=FT2KM;
maxalt*=FT2KM;
}
/* Get altitude */
if (coords_from_file && !arg_err && (alt < minalt || alt > maxalt))
{printf("\nError: unrecognized altitude %s in coordinate file line %1d\n\n",args[4],iline); arg_err = 1;}
while ((alt<minalt)||(alt>maxalt))
{
if (igdgc==2) printf("\n\nEnter geocentric altitude in km (%.2f to %.2f): ", minalt, maxalt);
if (igdgc==1 && units==1) printf("\n\nEnter geodetic altitude above mean sea level in km (%.2f to %.2f): ", minalt, maxalt);
if (igdgc==1 && units==2) printf("\n\nEnter geodetic altitude above mean sea level in meters (%.2f to %.2f): ", minalt, maxalt);
if (igdgc==1 && units==3) printf("\n\nEnter geodetic altitude above mean sea level in feet (%.2f to %.2f): ", minalt, maxalt);
safegets(inbuff, MAXREAD);
alt=atof(inbuff);
}
/* Convert altitude to km */
if (units==2)
{
alt *= 0.001;
}
else if (units==3)
{
alt /= FT2KM;
}
/* Get lat/long prefs */
if (coords_from_file && !arg_err && (decdeg != 1 && decdeg != 2))
{printf("\nError: unrecognized lat %s or lon %s in coordinate file line %1d\n\n",args[5],args[6],iline); arg_err = 1;}
while ((decdeg!=1)&&(decdeg!=2))
{
printf("\n\nHow would you like to enter the latitude and longitude?:");
printf("\n 1) In decimal degrees.");
printf("\n 2) In degrees, minutes, and seconds.\n");
printf("\n ==> ");
safegets(inbuff, MAXREAD);
decdeg=atoi(inbuff);
}
/* Get lat/lon */
if (decdeg==1)
{
if (coords_from_file && !arg_err && (latitude < -90 || latitude > 90))
{printf("\nError: unrecognized latitude %s in coordinate file line %1d\n\n",args[6],iline); arg_err = 1;}
while ((latitude<-90)||(latitude>90))
{
printf("\n\nEnter the decimal latitude (-90 to 90) (- for Southern hemisphere).\n");
safegets(inbuff, MAXREAD);
latitude=atof(inbuff);
}
if (coords_from_file && !arg_err && (longitude < -180 || longitude > 180))
{printf("\nError: unrecognized longitude %s in coordinate file line %1d\n\n",args[6],iline); arg_err = 1;}
while ((longitude<-180)||(longitude>180))
{
printf("\n\nEnter the decimal longitude (-180 to 180) (- for Western hemisphere).\n");
safegets(inbuff, MAXREAD);
longitude=atof(inbuff);
}
} /* if (decdeg==1) */
else
{
latitude=degrees_to_decimal(ilat_deg,ilat_min,ilat_sec);
longitude=degrees_to_decimal(ilon_deg,ilon_min,ilon_sec);
if (coords_from_file && !arg_err && (latitude < -90 || latitude > 90))
{printf("\nError: unrecognized latitude %s in coordinate file line %1d\n\n",args[6],iline); arg_err = 1;}
while ((latitude<-90)||(latitude>90))
{
ilat_deg=ilat_min=ilat_sec=200;
printf("\n\nEnter the decimal latitude (-90 to 90) (- for Southern hemisphere).\n");
while ((ilat_deg<-90)||(ilat_deg>90))
{
printf("\nDegrees (-90 to 90): ");
safegets(inbuff, MAXREAD);
ilat_deg=atoi(inbuff);
}
while ((ilat_min<-59)||(ilat_min>59))
{
printf("\nMinutes (-59 to 59): ");
safegets(inbuff, MAXREAD);
ilat_min=atoi(inbuff);
}
while ((ilat_sec<-59)||(ilat_sec>59))
{
printf("\nSeconds (-59 to 59): ");
safegets(inbuff, MAXREAD);
ilat_sec=atoi(inbuff);
}
latitude=degrees_to_decimal(ilat_deg,ilat_min,ilat_sec);
if ((latitude<-90)||(latitude>90))
{
printf("\nThe latitude %3.2f is out of range", latitude);
}
} /* while ((latitude<-90)||(latitude>90)) */
if (coords_from_file && !arg_err && (longitude < -180 || longitude > 180))
{printf("\nError: unrecognized longitude %s in coordinate file line %1d\n\n",args[6],iline); arg_err = 1;}
while ((longitude<-180)||(longitude>180))
{
ilon_deg=ilon_min=ilon_sec=200;
printf("\n\nEnter the decimal longitude (-180 to 180) (- for Western hemisphere).\n");
while ((ilon_deg<-180)||(ilon_deg>180))
{
printf("\nDegrees (-180 to 180): ");
safegets(inbuff, MAXREAD);
ilon_deg=atoi(inbuff);
}
while ((ilon_min<-59)||(ilon_min>59))
{
printf("\nMinutes (0 to 59): ");
safegets(inbuff, MAXREAD);
ilon_min=atoi(inbuff);
}
while ((ilon_sec<-59)||(ilon_sec>59))
{
printf("\nSeconds (0 to 59): ");
safegets(inbuff, MAXREAD);
ilon_sec=atoi(inbuff);
}
longitude=degrees_to_decimal(ilon_deg,ilon_min,ilon_sec);
if ((longitude<-180)||(longitude>180))
{
printf("\nThe longitude %3.2f is out of range", longitude);
}
} /* while ((longitude<-180)||(longitude>180)) */
} /* if (decdeg != 1) */
/** This will compute everything needed for 1 point in time. **/
if (max2[modelI] == 0)
{
getshc(mdfile, 1, irec_pos[modelI], max1[modelI], 1);
getshc(mdfile, 1, irec_pos[modelI+1], max1[modelI+1], 2);
nmax = interpsh(sdate, yrmin[modelI], max1[modelI],
yrmin[modelI+1], max1[modelI+1], 3);
nmax = interpsh(sdate+1, yrmin[modelI] , max1[modelI],
yrmin[modelI+1], max1[modelI+1],4);
}
else
{
getshc(mdfile, 1, irec_pos[modelI], max1[modelI], 1);
getshc(mdfile, 0, irec_pos[modelI], max2[modelI], 2);
nmax = extrapsh(sdate, epoch[modelI], max1[modelI], max2[modelI], 3);
nmax = extrapsh(sdate+1, epoch[modelI], max1[modelI], max2[modelI], 4);
}
/* Do the first calculations */
shval3(igdgc, latitude, longitude, alt, nmax, 3,
IEXT, EXT_COEFF1, EXT_COEFF2, EXT_COEFF3);
dihf(3);
shval3(igdgc, latitude, longitude, alt, nmax, 4,
IEXT, EXT_COEFF1, EXT_COEFF2, EXT_COEFF3);
dihf(4);
ddot = ((dtemp - d)*RAD2DEG);
if (ddot > 180.0) ddot -= 360.0;
if (ddot <= -180.0) ddot += 360.0;
ddot *= 60.0;
idot = ((itemp - i)*RAD2DEG)*60;
d = d*(RAD2DEG); i = i*(RAD2DEG);
hdot = htemp - h; xdot = xtemp - x;
ydot = ytemp - y; zdot = ztemp - z;
fdot = ftemp - f;
/* deal with geographic and magnetic poles */
if (h < 100.0) /* at magnetic poles */
{
d = NaN;
ddot = NaN;
/* while rest is ok */
}
if (h < 1000.0)
{
warn_H = 0;
warn_H_strong = 1;
if (h<warn_H_strong_val) warn_H_strong_val = h;
}
else if (h < 5000.0 && !warn_H_strong)
{
warn_H = 1;
if (h<warn_H_val) warn_H_val = h;
}
if (90.0-fabs(latitude) <= 0.001) /* at geographic poles */
{
x = NaN;
y = NaN;
d = NaN;
xdot = NaN;
ydot = NaN;
ddot = NaN;
warn_P = 1;
warn_H = 0;
warn_H_strong = 0;
/* while rest is ok */
}
/** Above will compute everything for 1 point in time. **/
/* Output the final results. */
if (coords_from_file)
{
print_result_file(outfile, d, i, h, x, y, z, f,ddot,idot,hdot,xdot,ydot,zdot,fdot);
}
else
{
printf("\n\n\n Model: %s \n", model[modelI]);
if (decdeg==1)
{
printf(" Latitude: %4.2f deg\n", latitude);
printf(" Longitude: %4.2f deg\n", longitude);
}
else
{
printf(" Latitude: %d deg, %d min, %d sec\n",
ilat_deg,ilat_min, ilat_sec);
printf(" Longitude: %d deg, %d min, %d sec\n",
ilon_deg, ilon_min, ilon_sec);
}
printf(" Altitude: ");
if (units==1)
printf("%.2f km\n", alt);
else if (units==2)
printf("%.2f meters\n", alt*1000.0);
else
printf("%.2f ft\n", (alt*FT2KM));
if (range==1)
{
printf(" Date of Interest: ");
if (decyears==1)
printf(" %4.2f\n\n", sdate);
else
printf("%d-%d-%d (yyyy-mm-dd)\n\n", isyear, ismonth, isday);
print_header();
print_result(sdate,d, i, h, x, y, z, f);
print_long_dashed_line();
print_header_sv();
print_result_sv(sdate,ddot,idot,hdot,xdot,ydot,zdot,fdot);
print_dashed_line();
} /* if range == 1 */
else
{
printf(" Range of Interest: ");
if (decyears==1)
printf("%4.2f to %4.2f, step %4.2f\n\n",sdate, edate, step);
else
printf("%d-%d-%d to %d-%d-%d (yyyy-mm-dd), step %4.2f (years)\n\n", isyear, ismonth, isday, ieyear, iemonth, ieday, step);
print_header();
print_result(sdate,d, i, h, x, y, z, f);
for(syr=sdate+step;syr<(edate+step);syr+=step)
{
if ((syr>edate)&&(edate!=(syr-step)))
{
syr=edate;
print_long_dashed_line();
}
/* Do the calculations */
for (counter=0;counter<step;counter++)
{
if (max2[modelI] == 0){ /*If not last element in array */
if (syr>yrmin[modelI+1]){ /* And past model boundary */
modelI++; /* Get next model */
}
}
} /* for counter */
if (max2[modelI] == 0) /*If still not last element */
{
getshc(mdfile, 1, irec_pos[modelI], max1[modelI], 1);
getshc(mdfile, 1, irec_pos[modelI+1], max1[modelI+1], 2);
nmax = interpsh(syr, yrmin[modelI], max1[modelI],
yrmin[modelI+1], max1[modelI+1], 3);
nmax = interpsh(syr+1, yrmin[modelI] , max1[modelI],
yrmin[modelI+1], max1[modelI+1],4);
}
else
{
getshc(mdfile, 1, irec_pos[modelI], max1[modelI], 1);
getshc(mdfile, 0, irec_pos[modelI], max2[modelI], 2);
nmax = extrapsh(syr, epoch[modelI], max1[modelI],
max2[modelI], 3);
nmax = extrapsh(syr+1, epoch[modelI], max1[modelI],
max2[modelI], 4);
}
shval3(igdgc, latitude, longitude, alt, nmax, 3,
IEXT, EXT_COEFF1, EXT_COEFF2, EXT_COEFF3);
dihf(3);
shval3(igdgc, latitude, longitude, alt, nmax, 4,
IEXT, EXT_COEFF1, EXT_COEFF2, EXT_COEFF3);
dihf(4);
ddot = ((dtemp - d)*RAD2DEG);
if (ddot > 180.0) ddot -= 360.0;
if (ddot <= -180.0) ddot += 360.0;
ddot *= 60.0;
idot = ((itemp - i)*RAD2DEG)*60.;
d = d*(RAD2DEG); i = i*(RAD2DEG);
hdot = htemp - h; xdot = xtemp - x;
ydot = ytemp - y; zdot = ztemp - z;
fdot = ftemp - f;
/* deal with geographic and magnetic poles */
if (h < 100.0) /* at magnetic poles */
{
d = NaN;
ddot = NaN;
/* while rest is ok */
}
if (90.0-fabs(latitude) <= 0.001) /* at geographic poles */
{
x = NaN;
y = NaN;
d = NaN;
xdot = NaN;
ydot = NaN;
ddot = NaN;
warn_P = 1;
warn_H = 0;
warn_H_strong = 0;
/* while rest is ok */
}
print_result(syr, d, i, h, x, y, z, f);
} /* for syr */
print_long_dashed_line();
print_header_sv();
print_result_sv(edate,ddot,idot,hdot,xdot,ydot,zdot,fdot);
print_dashed_line();
} /* if range > 1 */
if (warn_H)
{
printf("\nWarning: The horizontal field strength at this location is only %6.1f nT\n",warn_H_val);
printf(" Compass readings have large uncertainties in areas where H is\n");
printf(" smaller than 5000 nT\n\n");
}
if (warn_H_strong)
{
printf("\nWarning: The horizontal field strength at this location is only %6.1f nT\n",warn_H_strong_val);
printf(" Compass readings have VERY LARGE uncertainties in areas where H is\n");
printf(" smaller than 1000 nT\n\n");
}
if (warn_P)
{
printf("\nWarning: Location is at geographic pole where X, Y, and declination are not computed\n\n");
}
} /* if not coords_from_file */
if (coords_from_file)
again = !feof(coordfile) && !arg_err;
else
{
if (argc==7) again = 0; /* run command line only once */
else
do
{
printf("\n Enter");
printf("\n 0) to quit.");
printf("\n 1) to select a new model input file.");
printf("\n 2) to compute for a new point using same data file.");
printf("\n\n ==> ");
safegets(inbuff, MAXREAD);
again=atoi(inbuff);
if (again == 1) { need_to_read_model = 1; stream = NULL;}
if (again == 2) { need_to_read_model = 0; again = 1;}
}
while (again != 0 && again != 1);
} /* if not coords_from_file */
if (again == 1)
{
/* Reset defaults to catch on all while loops */
igdgc=decyears=units=decdeg=-1;
ismonth=isday=isyear=sdate=edate=range=step=-1;
latitude=ilat_deg=ilat_min=ilat_sec=200;
longitude=ilon_deg=ilon_min=ilon_sec=200;
alt=-9999999;
argc = 1;
}
} /* while (again == 1) */
reached_EOF:
if (coords_from_file) printf("\n Processed %1d lines\n\n",iline);
if (coords_from_file && !feof(coordfile) && arg_err) printf("Terminated prematurely due to argument error in coordinate file\n\n");
if (argc<7)
{
printf("\nThe same result could have been generated with the following command:\n");
printf("%s %s ", args[0], mdfile);
if (range == 1)
{
if (decyears==1)
{
printf("%4.2f ", sdate);
}
else
{
printf("%d,%d,%d ", isyear,ismonth,isday);
}
}
else
{
if (decyears==1)
{
printf("%4.2f-%4.2f-%4.2f ", sdate,edate,step);
}
else
{
printf("%d,%d,%d-%d,%d,%d-%4.2f ", isyear,ismonth,isday,ieyear,iemonth,ieday,step);
}
}
if (igdgc==1)
{
printf("D ");
}
else
{
printf("C ");
}
if (units==1)
{
printf("K%.2f ", alt);
}
else if (units==2)
{
printf("M%.2f ", alt*1000.0);
}
else
{
printf("F%.2f ", (alt*FT2KM));
}
if (decdeg==1)
{
printf("%3.2f %4.2f\n\n", latitude, longitude);
}
else
{
printf("%d,%d,%d %d,%d,%d\n\n", ilat_deg, ilat_min, ilat_sec,
ilon_deg, ilon_min, ilon_sec);
}
}
return 0;
}
void print_dashed_line(void)
{
printf(" -------------------------------------------------------------------------------\n");
return;
}
void print_long_dashed_line(void)
{
printf(" - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
return;
}
void print_header(void)
{
print_dashed_line();
printf(" Date D I H X Y Z F\n");
printf(" (yr) (deg min) (deg min) (nT) (nT) (nT) (nT) (nT)\n");
return;
}
void print_result(double date, double d, double i, double h, double x, double y, double z, double f)
{
int ddeg,ideg;
double dmin,imin;
/* Change d and i to deg and min */
ddeg=(int)d;
dmin=(d-(double)ddeg)*60;
if (d > 0 && dmin >= 59.5)
{
dmin -= 60.0;
ddeg++;
}
if (d < 0 && dmin <= -59.5)
{
dmin += 60.0;
ddeg--;
}
if (ddeg!=0) dmin=fabs(dmin);
ideg=(int)i;
imin=(i-(double)ideg)*60;
if (i > 0 && imin >= 59.5)
{
imin -= 60.0;
ideg++;
}
if (i < 0 && imin <= -59.5)
{
imin += 60.0;
ideg--;
}
if (ideg!=0) imin=fabs(imin);
if (my_isnan(d))
{
if (my_isnan(x))
printf(" %4.2f NaN %4dd %3.0fm %8.1f NaN NaN %8.1f %8.1f\n",date,ideg,imin,h,z,f);
else
printf(" %4.2f NaN %4dd %3.0fm %8.1f %8.1f %8.1f %8.1f %8.1f\n",date,ideg,imin,h,x,y,z,f);
}
else
printf(" %4.2f %4dd %3.0fm %4dd %3.0fm %8.1f %8.1f %8.1f %8.1f %8.1f\n",date,ddeg,dmin,ideg,imin,h,x,y,z,f);
return;
} /* print_result */
void print_header_sv(void)
{
printf(" Date dD dI dH dX dY dZ dF\n");
printf(" (yr) (min/yr) (min/yr) (nT/yr) (nT/yr) (nT/yr) (nT/yr) (nT/yr)\n");
} /* print_header_sv */
void print_result_sv(double date, double ddot, double idot, double hdot, double xdot, double ydot, double zdot, double fdot)
{
if (my_isnan(ddot))
{
if (my_isnan(xdot))
printf(" %4.2f NaN %7.1f %8.1f NaN NaN %8.1f %8.1f\n",date,idot,hdot,zdot,fdot);
else
printf(" %4.2f NaN %7.1f %8.1f %8.1f %8.1f %8.1f %8.1f\n",date,idot,hdot,xdot,ydot,zdot,fdot);
}
else
printf(" %4.2f %7.1f %7.1f %8.1f %8.1f %8.1f %8.1f %8.1f\n",date,ddot,idot,hdot,xdot,ydot,zdot,fdot);
return;
} /* print_result_sv */
void print_result_file(FILE *outf, double d, double i, double h, double x, double y, double z, double f,
double ddot, double idot, double hdot, double xdot, double ydot, double zdot, double fdot)
{
int ddeg,ideg;
double dmin,imin;
/* Change d and i to deg and min */
ddeg=(int)d;
dmin=(d-(double)ddeg)*60;
if (ddeg!=0) dmin=fabs(dmin);
ideg=(int)i;
imin=(i-(double)ideg)*60;
if (ideg!=0) imin=fabs(imin);
if (my_isnan(d))
{
if (my_isnan(x))
fprintf(outf," NaN %4dd %2.0fm %8.1f NaN NaN %8.1f %8.1f",ideg,imin,h,z,f);
else
fprintf(outf," NaN %4dd %2.0fm %8.1f %8.1f %8.1f %8.1f %8.1f",ideg,imin,h,x,y,z,f);
}
else
fprintf(outf," %4dd %2.0fm %4dd %2.0fm %8.1f %8.1f %8.1f %8.1f %8.1f",ddeg,dmin,ideg,imin,h,x,y,z,f);
if (my_isnan(ddot))
{
if (my_isnan(xdot))
fprintf(outf," NaN %7.1f %8.1f NaN NaN %8.1f %8.1f\n",idot,hdot,zdot,fdot);
else
fprintf(outf," NaN %7.1f %8.1f %8.1f %8.1f %8.1f %8.1f\n",idot,hdot,xdot,ydot,zdot,fdot);
}
else
fprintf(outf," %7.1f %7.1f %8.1f %8.1f %8.1f %8.1f %8.1f\n",ddot,idot,hdot,xdot,ydot,zdot,fdot);
return;
} /* print_result_file */
/****************************************************************************/
/* */
/* Subroutine safegets */
/* */
/****************************************************************************/
/* */
/* Gets characters from stdin untill it has reached n characters or \n, */
/* whichever comes first. \n is converted to \0. */
/* */
/* Input: n - Integer number of chars */
/* *buffer - Character array ptr which can contain n+1 characters */
/* */
/* Output: size - integer size of sting in buffer */
/* */
/* Note: All strings will be null terminated. */
/* */
/* By: David Owens */
/* dio@ngdc.noaa.gov */
/****************************************************************************/
int safegets(char *buffer,int n){
char *ptr; /** ptr used for finding '\n' **/
buffer[0] = '\0';
fgets(buffer,n,stdin); /** Get n chars **/
buffer[n+1]='\0'; /** Set last char to null **/
ptr=strchr(buffer,'\n'); /** If string contains '\n' **/
if (ptr!=NULL){ /** If string contains '\n' **/
ptr[0]='\0'; /** Change char to '\0' **/
if (buffer[0] == '\0') printf("\n ... no entry ...\n");
}
return strlen(buffer); /** Return the length **/
}
/****************************************************************************/
/* */
/* Subroutine degrees_to_decimal */
/* */
/****************************************************************************/
/* */
/* Converts degrees,minutes, seconds to decimal degrees. */
/* */
/* Input: */
/* degrees - Integer degrees */
/* minutes - Integer minutes */
/* seconds - Integer seconds */
/* */
/* Output: */
/* decimal - degrees in decimal degrees */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 12, 1988 */
/* */
/****************************************************************************/
double degrees_to_decimal(int degrees,int minutes,int seconds)
{
double deg;
double min;
double sec;
double decimal;
deg = degrees;
min = minutes/60.0;
sec = seconds/3600.0;
decimal = fabs(sec) + fabs(min) + fabs(deg);
if (deg < 0) {
decimal = -decimal;
} else if (deg == 0){
if (min < 0){
decimal = -decimal;
} else if (min == 0){
if (sec<0){
decimal = -decimal;
}
}
}
return(decimal);
}
/****************************************************************************/
/* */
/* Subroutine julday */
/* */
/****************************************************************************/
/* */
/* Computes the decimal day of year from month, day, year. */
/* Supplied by Daniel Bergstrom */
/* */
/* References: */
/* */
/* 1. Nachum Dershowitz and Edward M. Reingold, Calendrical Calculations, */
/* Cambridge University Press, 3rd edition, ISBN 978-0-521-88540-9. */
/* */
/* 2. Claus Tøndering, Frequently Asked Questions about Calendars, */
/* Version 2.9, http://www.tondering.dk/claus/calendar.html */
/* */
/****************************************************************************/
double julday(month, day, year)
int month;
int day;
int year;
{
int days[12] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
int leap_year = (((year % 4) == 0) &&
(((year % 100) != 0) || ((year % 400) == 0)));
double day_in_year = (days[month - 1] + day + (month > 2 ? leap_year : 0));
return ((double)year + (day_in_year / (365.0 + leap_year)));
}
/****************************************************************************/
/* */
/* Subroutine getshc */
/* */
/****************************************************************************/
/* */
/* Reads spherical harmonic coefficients from the specified */
/* model into an array. */
/* */
/* Input: */
/* stream - Logical unit number */
/* iflag - Flag for SV equal to ) or not equal to 0 */
/* for designated read statements */
/* strec - Starting record number to read from model */
/* nmax_of_gh - Maximum degree and order of model */
/* */
/* Output: */
/* gh1 or 2 - Schmidt quasi-normal internal spherical */
/* harmonic coefficients */
/* */
/* FORTRAN */
/* Bill Flanagan */
/* NOAA CORPS, DESDIS, NGDC, 325 Broadway, Boulder CO. 80301 */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 15, 1988 */
/* */
/****************************************************************************/
int getshc(file, iflag, strec, nmax_of_gh, gh)
char file[PATH];
int iflag;
long int strec;
int nmax_of_gh;
int gh;
{
char inbuff[MAXINBUFF];
char irat[9];
int ii,m,n,mm,nn;
int ios;
int line_num;
double g,hh;
double trash;
stream = fopen(file, "rt");
if (stream == NULL)
{
printf("\nError on opening file %s", file);
}
else
{
ii = 0;
ios = 0;
fseek(stream,strec,SEEK_SET);
for ( nn = 1; nn <= nmax_of_gh; ++nn)
{
for (mm = 0; mm <= nn; ++mm)
{
if (iflag == 1)
{
fgets(inbuff, MAXREAD, stream);
sscanf(inbuff, "%d%d%lg%lg%lg%lg%s%d",
&n, &m, &g, &hh, &trash, &trash, irat, &line_num);
}
else
{
fgets(inbuff, MAXREAD, stream);
sscanf(inbuff, "%d%d%lg%lg%lg%lg%s%d",
&n, &m, &trash, &trash, &g, &hh, irat, &line_num);
}
if ((nn != n) || (mm != m))
{
ios = -2;
fclose(stream);
return(ios);
}
ii = ii + 1;
switch(gh)
{
case 1: gh1[ii] = g;
break;
case 2: gh2[ii] = g;
break;
default: printf("\nError in subroutine getshc");
break;
}
if (m != 0)
{
ii = ii+ 1;
switch(gh)
{
case 1: gh1[ii] = hh;
break;
case 2: gh2[ii] = hh;
break;
default: printf("\nError in subroutine getshc");
break;
}
}
}
}
}
fclose(stream);
return(ios);
}
/****************************************************************************/
/* */
/* Subroutine extrapsh */
/* */
/****************************************************************************/
/* */
/* Extrapolates linearly a spherical harmonic model with a */
/* rate-of-change model. */
/* */
/* Input: */
/* date - date of resulting model (in decimal year) */
/* dte1 - date of base model */
/* nmax1 - maximum degree and order of base model */
/* gh1 - Schmidt quasi-normal internal spherical */
/* harmonic coefficients of base model */
/* nmax2 - maximum degree and order of rate-of-change model */
/* gh2 - Schmidt quasi-normal internal spherical */
/* harmonic coefficients of rate-of-change model */
/* */
/* Output: */
/* gha or b - Schmidt quasi-normal internal spherical */
/* harmonic coefficients */
/* nmax - maximum degree and order of resulting model */
/* */
/* FORTRAN */
/* A. Zunde */
/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 16, 1988 */
/* */
/****************************************************************************/
int extrapsh(date, dte1, nmax1, nmax2, gh)
double date;
double dte1;
int nmax1;
int nmax2;
int gh;
{
int nmax;
int k, l;
int ii;
double factor;
factor = date - dte1;
if (nmax1 == nmax2)
{
k = nmax1 * (nmax1 + 2);
nmax = nmax1;
}
else
{
if (nmax1 > nmax2)
{
k = nmax2 * (nmax2 + 2);
l = nmax1 * (nmax1 + 2);
switch(gh)
{
case 3: for ( ii = k + 1; ii <= l; ++ii)
{
gha[ii] = gh1[ii];
}
break;
case 4: for ( ii = k + 1; ii <= l; ++ii)
{
ghb[ii] = gh1[ii];
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
nmax = nmax1;
}
else
{
k = nmax1 * (nmax1 + 2);
l = nmax2 * (nmax2 + 2);
switch(gh)
{
case 3: for ( ii = k + 1; ii <= l; ++ii)
{
gha[ii] = factor * gh2[ii];
}
break;
case 4: for ( ii = k + 1; ii <= l; ++ii)
{
ghb[ii] = factor * gh2[ii];
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
nmax = nmax2;
}
}
switch(gh)
{
case 3: for ( ii = 1; ii <= k; ++ii)
{
gha[ii] = gh1[ii] + factor * gh2[ii];
}
break;
case 4: for ( ii = 1; ii <= k; ++ii)
{
ghb[ii] = gh1[ii] + factor * gh2[ii];
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
return(nmax);
}
/****************************************************************************/
/* */
/* Subroutine interpsh */
/* */
/****************************************************************************/
/* */
/* Interpolates linearly, in time, between two spherical harmonic */
/* models. */
/* */
/* Input: */
/* date - date of resulting model (in decimal year) */
/* dte1 - date of earlier model */
/* nmax1 - maximum degree and order of earlier model */
/* gh1 - Schmidt quasi-normal internal spherical */
/* harmonic coefficients of earlier model */
/* dte2 - date of later model */
/* nmax2 - maximum degree and order of later model */
/* gh2 - Schmidt quasi-normal internal spherical */
/* harmonic coefficients of internal model */
/* */
/* Output: */
/* gha or b - coefficients of resulting model */
/* nmax - maximum degree and order of resulting model */
/* */
/* FORTRAN */
/* A. Zunde */
/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 17, 1988 */
/* */
/****************************************************************************/
int interpsh(date, dte1, nmax1, dte2, nmax2, gh)
double date;
double dte1;
int nmax1;
double dte2;
int nmax2;
int gh;
{
int nmax;
int k, l;
int ii;
double factor;
factor = (date - dte1) / (dte2 - dte1);
if (nmax1 == nmax2)
{
k = nmax1 * (nmax1 + 2);
nmax = nmax1;
}
else
{
if (nmax1 > nmax2)
{
k = nmax2 * (nmax2 + 2);
l = nmax1 * (nmax1 + 2);
switch(gh)
{
case 3: for ( ii = k + 1; ii <= l; ++ii)
{
gha[ii] = gh1[ii] + factor * (-gh1[ii]);
}
break;
case 4: for ( ii = k + 1; ii <= l; ++ii)
{
ghb[ii] = gh1[ii] + factor * (-gh1[ii]);
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
nmax = nmax1;
}
else
{
k = nmax1 * (nmax1 + 2);
l = nmax2 * (nmax2 + 2);
switch(gh)
{
case 3: for ( ii = k + 1; ii <= l; ++ii)
{
gha[ii] = factor * gh2[ii];
}
break;
case 4: for ( ii = k + 1; ii <= l; ++ii)
{
ghb[ii] = factor * gh2[ii];
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
nmax = nmax2;
}
}
switch(gh)
{
case 3: for ( ii = 1; ii <= k; ++ii)
{
gha[ii] = gh1[ii] + factor * (gh2[ii] - gh1[ii]);
}
break;
case 4: for ( ii = 1; ii <= k; ++ii)
{
ghb[ii] = gh1[ii] + factor * (gh2[ii] - gh1[ii]);
}
break;
default: printf("\nError in subroutine extrapsh");
break;
}
return(nmax);
}
/****************************************************************************/
/* */
/* Subroutine shval3 */
/* */
/****************************************************************************/
/* */
/* Calculates field components from spherical harmonic (sh) */
/* models. */
/* */
/* Input: */
/* igdgc - indicates coordinate system used; set equal */
/* to 1 if geodetic, 2 if geocentric */
/* latitude - north latitude, in degrees */
/* longitude - east longitude, in degrees */
/* elev - WGS84 altitude above ellipsoid (igdgc=1), or */
/* radial distance from earth's center (igdgc=2) */
/* a2,b2 - squares of semi-major and semi-minor axes of */
/* the reference spheroid used for transforming */
/* between geodetic and geocentric coordinates */
/* or components */
/* nmax - maximum degree and order of coefficients */
/* iext - external coefficients flag (=0 if none) */
/* ext1,2,3 - the three 1st-degree external coefficients */
/* (not used if iext = 0) */
/* */
/* Output: */
/* x - northward component */
/* y - eastward component */
/* z - vertically-downward component */
/* */
/* based on subroutine 'igrf' by D. R. Barraclough and S. R. C. Malin, */
/* report no. 71/1, institute of geological sciences, U.K. */
/* */
/* FORTRAN */
/* Norman W. Peddie */
/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 17, 1988 */
/* */
/****************************************************************************/
int shval3(igdgc, flat, flon, elev, nmax, gh, iext, ext1, ext2, ext3)
int igdgc;
double flat;
double flon;
double elev;
int nmax;
int gh;
int iext;
double ext1;
double ext2;
double ext3;
{
double earths_radius = 6371.2;
double dtr = 0.01745329;
double slat;
double clat;
double ratio;
double aa, bb, cc, dd;
double sd;
double cd;
double r;
double a2;
double b2;
double rr;
double fm,fn;
double sl[14];
double cl[14];
double p[119];
double q[119];
int ii,j,k,l,m,n;
int npq;
int ios;
double argument;
double power;
a2 = 40680631.59; /* WGS84 */
b2 = 40408299.98; /* WGS84 */
ios = 0;
r = elev;
argument = flat * dtr;
slat = sin( argument );
if ((90.0 - flat) < 0.001)
{
aa = 89.999; /* 300 ft. from North pole */
}
else
{
if ((90.0 + flat) < 0.001)
{
aa = -89.999; /* 300 ft. from South pole */
}
else
{
aa = flat;
}
}
argument = aa * dtr;
clat = cos( argument );
argument = flon * dtr;
sl[1] = sin( argument );
cl[1] = cos( argument );
switch(gh)
{
case 3: x = 0;
y = 0;
z = 0;
break;
case 4: xtemp = 0;
ytemp = 0;
ztemp = 0;
break;
default: printf("\nError in subroutine shval3");
break;
}
sd = 0.0;
cd = 1.0;
l = 1;
n = 0;
m = 1;
npq = (nmax * (nmax + 3)) / 2;
if (igdgc == 1)
{
aa = a2 * clat * clat;
bb = b2 * slat * slat;
cc = aa + bb;
argument = cc;
dd = sqrt( argument );
argument = elev * (elev + 2.0 * dd) + (a2 * aa + b2 * bb) / cc;
r = sqrt( argument );
cd = (elev + dd) / r;
sd = (a2 - b2) / dd * slat * clat / r;
aa = slat;
slat = slat * cd - clat * sd;
clat = clat * cd + aa * sd;
}
ratio = earths_radius / r;
argument = 3.0;
aa = sqrt( argument );
p[1] = 2.0 * slat;
p[2] = 2.0 * clat;
p[3] = 4.5 * slat * slat - 1.5;
p[4] = 3.0 * aa * clat * slat;
q[1] = -clat;
q[2] = slat;
q[3] = -3.0 * clat * slat;
q[4] = aa * (slat * slat - clat * clat);
for ( k = 1; k <= npq; ++k)
{
if (n < m)
{
m = 0;
n = n + 1;
argument = ratio;
power = n + 2;
rr = pow(argument,power);
fn = n;
}
fm = m;
if (k >= 5)
{
if (m == n)
{
argument = (1.0 - 0.5/fm);
aa = sqrt( argument );
j = k - n - 1;
p[k] = (1.0 + 1.0/fm) * aa * clat * p[j];
q[k] = aa * (clat * q[j] + slat/fm * p[j]);
sl[m] = sl[m-1] * cl[1] + cl[m-1] * sl[1];
cl[m] = cl[m-1] * cl[1] - sl[m-1] * sl[1];
}
else
{
argument = fn*fn - fm*fm;
aa = sqrt( argument );
argument = ((fn - 1.0)*(fn-1.0)) - (fm * fm);
bb = sqrt( argument )/aa;
cc = (2.0 * fn - 1.0)/aa;
ii = k - n;
j = k - 2 * n + 1;
p[k] = (fn + 1.0) * (cc * slat/fn * p[ii] - bb/(fn - 1.0) * p[j]);
q[k] = cc * (slat * q[ii] - clat/fn * p[ii]) - bb * q[j];
}
}
switch(gh)
{
case 3: aa = rr * gha[l];
break;
case 4: aa = rr * ghb[l];
break;
default: printf("\nError in subroutine shval3");
break;
}
if (m == 0)
{
switch(gh)
{
case 3: x = x + aa * q[k];
z = z - aa * p[k];
break;
case 4: xtemp = xtemp + aa * q[k];
ztemp = ztemp - aa * p[k];
break;
default: printf("\nError in subroutine shval3");
break;
}
l = l + 1;
}
else
{
switch(gh)
{
case 3: bb = rr * gha[l+1];
cc = aa * cl[m] + bb * sl[m];
x = x + cc * q[k];
z = z - cc * p[k];
if (clat > 0)
{
y = y + (aa * sl[m] - bb * cl[m]) *
fm * p[k]/((fn + 1.0) * clat);
}
else
{
y = y + (aa * sl[m] - bb * cl[m]) * q[k] * slat;
}
l = l + 2;
break;
case 4: bb = rr * ghb[l+1];
cc = aa * cl[m] + bb * sl[m];
xtemp = xtemp + cc * q[k];
ztemp = ztemp - cc * p[k];
if (clat > 0)
{
ytemp = ytemp + (aa * sl[m] - bb * cl[m]) *
fm * p[k]/((fn + 1.0) * clat);
}
else
{
ytemp = ytemp + (aa * sl[m] - bb * cl[m]) *
q[k] * slat;
}
l = l + 2;
break;
default: printf("\nError in subroutine shval3");
break;
}
}
m = m + 1;
}
if (iext != 0)
{
aa = ext2 * cl[1] + ext3 * sl[1];
switch(gh)
{
case 3: x = x - ext1 * clat + aa * slat;
y = y + ext2 * sl[1] - ext3 * cl[1];
z = z + ext1 * slat + aa * clat;
break;
case 4: xtemp = xtemp - ext1 * clat + aa * slat;
ytemp = ytemp + ext2 * sl[1] - ext3 * cl[1];
ztemp = ztemp + ext1 * slat + aa * clat;
break;
default: printf("\nError in subroutine shval3");
break;
}
}
switch(gh)
{
case 3: aa = x;
x = x * cd + z * sd;
z = z * cd - aa * sd;
break;
case 4: aa = xtemp;
xtemp = xtemp * cd + ztemp * sd;
ztemp = ztemp * cd - aa * sd;
break;
default: printf("\nError in subroutine shval3");
break;
}
return(ios);
}
/****************************************************************************/
/* */
/* Subroutine dihf */
/* */
/****************************************************************************/
/* */
/* Computes the geomagnetic d, i, h, and f from x, y, and z. */
/* */
/* Input: */
/* x - northward component */
/* y - eastward component */
/* z - vertically-downward component */
/* */
/* Output: */
/* d - declination */
/* i - inclination */
/* h - horizontal intensity */
/* f - total intensity */
/* */
/* FORTRAN */
/* A. Zunde */
/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */
/* */
/* C */
/* C. H. Shaffer */
/* Lockheed Missiles and Space Company, Sunnyvale CA */
/* August 22, 1988 */
/* */
/****************************************************************************/
int dihf (gh)
int gh;
{
int ios;
int j;
double sn;
double h2;
double hpx;
double argument, argument2;
ios = gh;
sn = 0.0001;
switch(gh)
{
case 3: for (j = 1; j <= 1; ++j)
{
h2 = x*x + y*y;
argument = h2;
h = sqrt(argument); /* calculate horizontal intensity */
argument = h2 + z*z;
f = sqrt(argument); /* calculate total intensity */
if (f < sn)
{
d = NaN; /* If d and i cannot be determined, */
i = NaN; /* set equal to NaN */
}
else
{
argument = z;
argument2 = h;
i = atan2(argument,argument2);
if (h < sn)
{
d = NaN;
}
else
{
hpx = h + x;
if (hpx < sn)
{
d = PI;
}
else
{
argument = y;
argument2 = hpx;
d = 2.0 * atan2(argument,argument2);
}
}
}
}
break;
case 4: for (j = 1; j <= 1; ++j)
{
h2 = xtemp*xtemp + ytemp*ytemp;
argument = h2;
htemp = sqrt(argument);
argument = h2 + ztemp*ztemp;
ftemp = sqrt(argument);
if (ftemp < sn)
{
dtemp = NaN; /* If d and i cannot be determined, */
itemp = NaN; /* set equal to 999.0 */
}
else
{
argument = ztemp;
argument2 = htemp;
itemp = atan2(argument,argument2);
if (htemp < sn)
{
dtemp = NaN;
}
else
{
hpx = htemp + xtemp;
if (hpx < sn)
{
dtemp = PI;
}
else
{
argument = ytemp;
argument2 = hpx;
dtemp = 2.0 * atan2(argument,argument2);
}
}
}
}
break;
default: printf("\nError in subroutine dihf");
break;
}
return(ios);
}
Reference in New Issue Copy Permalink