/* COMPILATION: cc -lm geonew.c -o geonew /****************************************************************************/ /* */ /* Disclaimer: This C version of the Geomagnetic Field */ /* Modeling software is being supplied to aid programmers in */ /* integrating spherical harmonic field modeling calculations */ /* into their code. It is being distributed unoffically. The */ /* National Geophysical Data Center does not support it as */ /* a data product or guarantee it's correctness. The limited */ /* testing done on this code seems to indicate that the results */ /* obtained from this program are compariable to those obtained */ /* */ /****************************************************************************/ /* */ /* In regards to the disclaimer, to the best of my knowlege and quick */ /* testing, this program, generates numaric output which is withing .1 */ /* degress of the current fortran version. However, it *is* a program */ /* and most likely contains bugs. */ /* dio 6-6-96 */ /* */ /****************************************************************************/ /* */ /* This is version 3.01 of the source code but still represents the */ /* geomag30 executable. */ /* dio 9-17-96 */ /* */ /****************************************************************************/ /* */ /* Bug fix- the range of dates error. There was a difference in the */ /* output between the last value on a range of dates and the individual */ /* value for the same date. Lets make this version 3.02 */ /* ljh 11-20-98 */ /* */ /****************************************************************************/ /* Program was modified so that it can accept year 2000 models */ /* This required that the number of blank spaces on model header */ /* records be decreased from four to three */ /* */ /* 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, */ /* elevation, 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 (measured from */ /* the surface of the earth), */ /* or geocentric (measured from the */ /* center of the earth). */ /* */ /* Elevation : Elevation above sea level in kilometers.*/ /* if geocentric coordinate preference is */ /* used then the elevation must be in the */ /* range of 6370.20 km - 6971.20 km as */ /* measured from the center of the earth. */ /* Enter elevation in kilometers in the */ /* form xxxx.xx. */ /* */ /* 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 #include #include /* standard C include statements */ #include #include #include #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 (float)0 #define EXT_COEFF2 (float)0 #define EXT_COEFF3 (float)0 float gh1[171]; float gh2[171]; float gha[171]; /* Geomag global variables */ float ghb[171]; float d=0,f=0,h=0,i=0; float dtemp,ftemp,htemp,itemp; float x=0,y=0,z=0; float 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 */ /* Lockheed Missiles and Space Company */ /* Sunnyvale, California */ /* (408) 756 - 5210 */ /* 29 July, 1988 */ /* */ /* Contact : Greg Fountain */ /* Lockheed Missiles and Space Company */ /* Dept. MO-30, Bldg. 581 */ /* 1111 Lockheed Way */ /* P.O. Box 3504 */ /* Sunnyvale, Calif. 94088-3504 */ /* */ /* Rewritten by : David Owens */ /* dio@ngdc.noaa.gov */ /* For Susan McClean */ /* */ /* Contact : Susan McLean */ /* sjm@ngdc.noaa.gov */ /* National Geophysical Data Center */ /* World Data Center-A for Solid Earth Geophysics */ /* NOAA, E/GC1, 325 Broadway, */ /* Boulder, CO 80303 */ /* */ /* Original */ /* FORTRAN */ /* Programmer : National Geophysical Data Center */ /* World Data Center-A for Solid Earth Geophysics */ /* NOAA, E/GC1, 325 Broadway, */ /* Boulder, CO 80303 */ /* */ /* Contact : National Geophysical Data Center */ /* World Data Center-A for Solid Earth Geophysics */ /* NOAA, E/GC1, 325 Broadway, */ /* Boulder, CO 80303 */ /* */ /****************************************************************************/ /* */ /* dio modifications include overhauling the interactive interface to */ /* support platform independence and improve fatal error dectection and */ /* prevention. A command line interface was added and the interactive */ /* interface was streamlined. The function safegets() was added and the */ /* function getshc's i/0 was modified. A option to specify a date range */ /* range and step was also added. */ /* */ /* dio 6-6-96 */ /* */ /****************************************************************************/ /* */ /* Some variables used in this program */ /* */ /* Name Type Usage */ /* ------------------------------------------------------------------------ */ /* */ /* a2,b2 Scalar Float Squares of semi-major and semi-minor */ /* axes of the reference spheroid used */ /* for transforming between geodetic or */ /* geocentric coordinates. */ /* */ /* minalt Float array of MAXMOD Minimum height of model. */ /* */ /* altmin Float Minimum height of selected model. */ /* */ /* altmax Float array of MAXMOD Maximum height of model. */ /* */ /* maxalt Float Maximum height of selected model. */ /* */ /* d Scalar Float Declination of the field from the */ /* geographic north (deg). */ /* */ /* sdate Scalar Float start date inputted */ /* */ /* ddot Scalar Float annual rate of change of decl. */ /* (deg/yr) */ /* */ /* elev Scalar Float elevation. */ /* */ /* epoch Float array of MAXMOD epoch of model. */ /* */ /* ext Scalar Float Three 1st-degree external coeff. */ /* */ /* latitude Scalar Float Latitude. */ /* */ /* longitude Scalar Float Longitude. */ /* */ /* gh1 Float array of 120 Schmidt quasi-normal internal */ /* spherical harmonic coeff. */ /* */ /* gh2 Float array of 120 Schmidt quasi-normal internal */ /* spherical harmonic coeff. */ /* */ /* gha Float array of 120 Coefficients of resulting model. */ /* */ /* ghb Float array of 120 Coefficients of rate of change model.*/ /* */ /* i Scalar Float Inclination (deg). */ /* */ /* idot Scalar Float Rate of change of i (deg/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 Float Min year of all models */ /* */ /* maxyr Float Max year of all models */ /* */ /* yrmax Float array of MAXMOD Max year of model. */ /* */ /* yrmin Float array of MAXMOD Min year of model. */ /* */ /****************************************************************************/ int main(int argv, char**argc) { #ifdef MAC ccommand(argv, argc); #endif /* Variable declaration */ /* Control variables */ int again = 1; int decyears = 3; int units = 4; int decdeg = 3; int range = -1; int counter = 0; int modelI; /* Which model (Index) */ 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; int ddeg, ideg; long irec_pos[MAXMOD]; char *mdfile; char inbuff[MAXINBUFF]; char model[MAXMOD][8]; char *fila; char *begin; char *rest; FILE *cano; float epoch[MAXMOD]; float yrmin[MAXMOD]; float yrmax[MAXMOD]; float minyr; float maxyr; float altmin[MAXMOD]; float altmax[MAXMOD]; float minalt; float maxalt; float elev=-999999; float sdate=-1; float step=-1; float syr; float edate=-1; float latitude=200; float longitude=200; float ddot; float fdot; float hdot; float idot; float xdot; float in_field; float out_field; float ydot; float zdot; float dmin, imin; /* Subroutines used */ float degrees_to_decimal(); float julday(); int interpsh(); int extrapsh(); int shval3(); int dihf(); int safegets(char *buffer,int n); /* Initializations. */ inbuff[MAXREAD+1]='\0'; /* Just to protect mem. */ inbuff[MAXINBUFF-1]='\0'; /* Just to protect mem. */ mdfile = "igrf"; stream=fopen(mdfile, "rt"); /* Obtain the desired model file and read the data */ 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%f%d%d%d%f%f%f%f", 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]maxyr){ maxyr=yrmax[modelI]; } } } } fclose(stream); /* Take in field data */ fila = "t.in"; cano = fopen(fila, "rt"); for (i=0; i<30; i=1) { if(fscanf(cano, "%d %d %d %f %f %f", &isyear, &ismonth, &isday, &latitude, &longitude, &in_field) == EOF) {break;} else { /* Get date */ sdate = julday(ismonth,isday,isyear); elev = 0.; igdgc = 1; /* Pick model */ for(modelI=0;sdate>yrmax[modelI];modelI++); if(sdate 0.0) && (i_month > 2)) { leap_year = 1; } else { leap_year = 0; } } day = aggregate_first_day_of_month[i_month] + i_day - 1 + leap_year; if (leap_year) { decimal_date = year + (day/366.0); /*In version 3.0 this was incorrect*/ } else { decimal_date = year + (day/365.0); /*In version 3.0 this was incorrect*/ } return(decimal_date); } /****************************************************************************/ /* */ /* 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; float g,hh; float 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, 3, stream); inbuff[3]='\0'; sscanf(inbuff, "%d", &m); fgets(inbuff, 3, stream); inbuff[3]='\0'; sscanf(inbuff, "%d", &n); fgets(inbuff, MAXREAD-4, stream); sscanf(inbuff, "%f%f%f%f%s%d", &g, &hh, &trash, &trash, irat, &line_num); } else { fgets(inbuff, 3, stream); inbuff[3]='\0'; sscanf(inbuff, "%d", &m); fgets(inbuff, 3, stream); inbuff[3]='\0'; sscanf(inbuff, "%d", &n); fgets(inbuff, MAXREAD-4, stream); sscanf(inbuff, "%f%f%f%f%s%d", &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) float date; float dte1; int nmax1; int nmax2; int gh; { int nmax; int k, l; int ii; float 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) float date; float dte1; int nmax1; float dte2; int nmax2; int gh; { int nmax; int k, l; int ii; float 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 - elevation above mean sea level (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; float flat; float flon; float elev; int nmax; int gh; int iext; float ext1; float ext2; float ext3; { float earths_radius = 6371.2; float dtr = 0.01745329; float slat; float clat; float ratio; float aa, bb, cc, dd; float sd; float cd; float r; float a2; float b2; float rr; float fm,fn; float sl[14]; float cl[14]; float p[119]; float q[119]; int ii,j,k,l,m,n; int npq; int ios; double arguement; double power; a2 = 40680925.0; b2 = 40408588.0; ios = 0; r = elev; arguement = flat * dtr; slat = sin( arguement ); 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; } } arguement = aa * dtr; clat = cos( arguement ); arguement = flon * dtr; sl[1] = sin( arguement ); cl[1] = cos( arguement ); 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; arguement = cc; dd = sqrt( arguement ); arguement = elev * (elev + 2.0 * dd) + (a2 * aa + b2 * bb) / cc; r = sqrt( arguement ); 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; arguement = 3.0; aa = sqrt( arguement ); 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; arguement = ratio; power = n + 2; rr = pow(arguement,power); fn = n; } fm = m; if (k >= 5) { if (m == n) { arguement = (1.0 - 0.5/fm); aa = sqrt( arguement ); 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 { arguement = fn*fn - fm*fm; aa = sqrt( arguement ); arguement = ((fn - 1.0)*(fn-1.0)) - (fm * fm); bb = sqrt( arguement )/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; float sn; float h2; float hpx; double arguement, arguement2; double rad, pi; ios = gh; sn = 0.0001; rad = 57.29577951; pi = 3.141592654; switch(gh) { case 3: for (j = 1; j <= 1; ++j) { h2 = x*x + y*y; arguement = h2; h = sqrt(arguement); /* calculate horizontal intensity */ arguement = h2 + z*z; f = sqrt(arguement); /* calculate total intensity */ if (f < sn) { d = 999.0 * rad; /* If d and i cannot be determined, */ i = 999.0 * rad; /* set equal to 999.0 */ } else { arguement = z; arguement2 = h; i = atan2(arguement,arguement2); if (h < sn) { d = 999.0 * rad; } else { hpx = h + x; if (hpx < sn) { d = pi; } else { arguement = y; arguement2 = hpx; d = 2.0 * atan2(arguement,arguement2); } } } } break; case 4: for (j = 1; j <= 1; ++j) { h2 = xtemp*xtemp + ytemp*ytemp; arguement = h2; htemp = sqrt(arguement); arguement = h2 + ztemp*ztemp; ftemp = sqrt(arguement); if (ftemp < sn) { dtemp = 999.0; /* If d and i cannot be determined, */ itemp = 999.0; /* set equal to 999.0 */ } else { arguement = ztemp; arguement2 = htemp; itemp = atan2(arguement,arguement2); if (htemp < sn) { dtemp = 999.0; } else { hpx = htemp + xtemp; if (hpx < sn) { dtemp = 180.0; } else { arguement = ytemp; arguement2 = hpx; dtemp = 2.0 * atan2(arguement,arguement2); } } } } break; default: printf("\nError in subroutine dihf"); break; } return(ios); }