C 8/92: /home/cedar/emery/geomag/psiwt3.f C Input file is filein. Output is flout2 with corrected psi values C in same format, and date.psi (flout1) in format for program convert. C Corrected geomagnetic coordinates are also found or corrected. C 8/92: Apex coordinates used (Gang Lu), plus only picks out selected C stations from master list with iget. C To compile: f77 -o psiwt3.ob psiwt3.f apex.f PROGRAM PSIWT3 character*90 line character*19 stat character*3 stat3,comp character*4 statnm character*80 flout1 character*80 filein character*20 period dimension iget(700),iskp(700) data iget/700*1/ integer nsta date = 1991.58 period = 'August 2-3, 1991' nsta = 103 nday = 2 minsp = 5 write(6,*) 'Enter file to read :' read(5,'(A80)') filein write(6,*) 'Enter file to output :' read(5,'(A80)') flout1 write(6,*) 'Enter period (ex. August 2-3, 1991):' read(5,'(A20)') period write(6,*) 'Enter number of days :' read(5,*) nday write(6,*) 'Enter time over which data is averaged (in minutes)' read(5,*) minsp write(6,*) 'filein,fileout,period = ',filein,flout1,period write(6,*) 'Enter date (ex 1991.85) : ' read(5,*) date write(6,*) 'Date = ',date open (1,file=filein,status='old') read(1,"(I3)") nsta C ********************** Will change ****************************** C Assume have full days numtm = nday*24*60/minsp C Check to see if # of iget=1 is the same as nsta niget = 0 niskp = 0 do 90 i=1,700 if (iget(i) .eq. 1) niget = niget + 1 if (iskp(i) .eq. 1) niskp = niskp + 1 90 continue write (6,"(1x,'nsta iget=1 niskp: ',3i4)") nsta,niget,niskp if (niget .ne. nsta) write (6,"(1x,'# of iget=1 is NOT' | ' the same as nsta -- WARNING')") C ***add a subroutine to catch overfloating *** C call fp_abort open (11,file=flout1,status='unknown') write (11,"(i3,' 0 0',i5,i3,' nsta minxy minz ntm minsp ', | a15)") nsta,numtm,minsp,period write (11,"(' # name cmp h a glat glon alat alon', | ' dec psi wt')") C Read 3 dummy lines do 10 k = 1,2 read (1,"(a90)") line 10 continue j = 1 k = 1 m = 1 100 READ(1,1000,END=150)stat3,stat,statnm,comp,glat,glon,bLAT,bLON,dec 1000 FORMAT(a3,1x,a19,a4,4x,a3,2x,2f7.2,1x,2f7.2,2x,f7.2) alt = 110 call apex(date,glat,glon,alt,1,a,alat,alon, | psi,bmag,xmag,ymag,zmag,clate,elon) ang = psi wt = 1.0 if (abs(alat).gt.50.0) wt = 0.5 if (abs(alat).gt.80.0) wt = 0.25 if (statnm.eq.'ALER') wt = 0.12 if (statnm.eq.'HEIS') wt = 0.25 1001 FORMAT(i3,1x,a3,1x,a19,a4,4x,a3,2x,2f7.2,1x,2f7.2,2x,2f7.2) k = k + 1 if (iget(k-1) .eq. 0) go to 100 WRITE (11,1002) j,stat3,statnm,comp,glat,glon,alat,alon,dec,ang,wt 1002 FORMAT(i3,1x,a3,1x,a4,1x,a3,4x,2F7.2,1x,2f7.2,2x,3f7.2) j = j + 1 if (iskp(k-1) .eq. 1) go to 100 m = m + 1 GO TO 100 150 STOP END C************************************************************************* C CG SP: (-74.3,127.5) in geog coords (colat=164.3 or 15.7 if like NH) C INPUTS: C CLAT = CORRECTED GEOMAGNETIC LATITUDE OF STATION C GLON = GEOGRAPHIC LONGITUDE OF STATION C OUTPUT: C ANG = ANGLE BETWEEN CG south AND GEOGR. south (PSI) C ang is positive if go clockwise (looking down!!!) between line C joining station and SP and line joining station and CG SP. C************************************************************************* C EQUATION: C SIN(-PSI) = SIN(CNPCOL)*(SIN(GLON-CNPLON))/SIN(CCOLAT) C WHERE: C PSI IS OUTPUT ANGLE (ANG) C CNPCOL IS GEOGRAPHIC COLATITUDE OF THE CORR. GEOMAG. south POLE C assuming it is in the NH. (i.e., 15.7 deg) C GLON IS GEOGRAPHIC EAST LONGITUDE OF STATION (INPUT) C CNPLON IS GEOGRAPHIC EAST LONGITUDE OF THE CORR. GEOMAG. south C pole + 180 deg to convert calc to NH. (i.e. 307.5) C CCOLAT IS CORR. GEOMAG. COLATITUDE OF STATION (FOUND FROM C CORR. GEOMAG. LATITUDE OF STATION (INPUT) ) C************************************************************************** SUBROUTINE shpsi (CLAT,GLON,ANG,j,stat3) PARAMETER ( DGTORD = 3.1415926535898/180.0 ) C CHANGE CG LAT INTO CG COLAT FOR STATION C AND CHANGE ALL ANGLES FROM DEGREES TO RADS write(6,"('shpsi: CLAT,GLON,j,stat3',2f8.2,i3,a3)") | CLAT,GLON,j,stat3 CCOLAT = (90.0 - abs(CLAT))*DGTORD CNPCOL = 15.7*DGTORD CNPLON = 307.5*DGTORD GGLON = GLON*DGTORD C CALCULATE NUMERATOR AND DENOMINATOR OF EQN OPP = SIN(CNPCOL) * SIN(GGLON-CNPLON) HYPOT = SIN(CCOLAT) C FIND SIN(-PSI) SPSI = OPP/HYPOT if (abs(SPSI) .ge. 1.) then write(6,"('stat3,j,CLAT,GLON,CCOLAT,CNPCOL,CNPLON,GGLON,' |,'OPP,HYPOT=',a3,i3,8f8.4)")stat3,j,CLAT,GLON,CCOLAT,CNPCOL, | CNPLON,GGLON,OPP,HYPOT write(6,"('SPSI=',f5.2)")SPSI end if ANG = - ASIN(SPSI) C CHANGE ANGLE BACK TO DEGREES ANG = ANG/DGTORD RETURN END c ************************************************************************* C************************************************************************* C SUBROUTINE CONVERT (CLAT,GLON,ANG) C THIS SUBROUTINE FINDS THE ANGLE BETWEEN CORRECTED GEOMAGNETIC C NORTH AND GEOGRAPHIC NORTH FOR A GIVEN LOCATION. THE ANGLE IS C POSITIVE WHEN THE CORRECTED GEOMAGNETIC POLE IS EASTWARD OF THE C GEOGRAPHIC POLE, AND NEGATIVE WHEN THE CG POLE IS WESTWARD. C WRITTEN BY J. MCKEE, NAT.L GEOPHYSICAL DATA CENTER, C 325 BROADWAY, BOULDER, CO 80303 NOVEMBER 1989 C************************************************************************* C INPUTS: C CLAT = CORRECTED GEOMAGNETIC LATITUDE OF STATION C GLON = GEOGRAPHIC LONGITUDE OF STATION C OUTPUT: C ANG = ANGLE BETWEEN CG NORTH AND GEOGR. NORTH (PSI) C ang is positive if go clockwise from line between station and NP C to line between station and CG NP. C************************************************************************* C EQUATION: C SIN(-PSI) = SIN(CNPCOL)*(SIN(GLON-CNPLON))/SIN(CCOLAT) C WHERE: C PSI IS OUTPUT ANGLE (ANG) C CNPCOL IS GEOGRAPHIC COLATITUDE OF THE CORR. GEOMAG. NORTH POLE C (9.38 DEGREES) C GLON IS GEOGRAPHIC EAST LONGITUDE OF STATION (INPUT) C CNPLON IS GEOGRAPHIC EAST LONGITUDE OF THE CORR. GEOMAG. NORTH C POLE (280.12 DEGREES) C CCOLAT IS CORR. GEOMAG. COLATITUDE OF STATION (FOUND FROM C CORR. GEOMAG. LATITUDE OF STATION (INPUT) ) C************************************************************************** SUBROUTINE CONVERT (CLAT,GLON,ANG) PARAMETER ( DGTORD = 3.1415926535898/180.0 ) C CHANGE CG LAT INTO CG COLAT FOR STATION C AND CHANGE ALL ANGLES FROM DEGREES TO RADS CCOLAT = (90.0 - CLAT)*DGTORD CNPCOL = 9.38*DGTORD CNPLON = 280.12*DGTORD GGLON = GLON*DGTORD C CALCULATE NUMERATOR AND DENOMINATOR OF EQN OPP = SIN(CNPCOL) * SIN(GGLON-CNPLON) HYPOT = SIN(CCOLAT) C FIND SIN(-PSI) SPSI = OPP/HYPOT IF (ABS(SPSI) .GE. 1.) WRITE (6,"('SPSI=',F5.2)") SPSI ANG = - ASIN(SPSI) C CHANGE ANGLE BACK TO DEGREES ANG = ANG/DGTORD IF (ANG .LT. -180.) ANG = 360. + ANG RETURN END c ************************************************************************* subroutine cnvcoord(inlat,inlong,height,order, & outlat,outlong,outr,mgflag,err) c c This routine uses a set of coefficients to perform a coordinate c conversion between geographic and magnetic (or visa versa) c coordinates. The conversion is done using a spherical c harmonic expansion. c parameter ( Lorder=25,pi=3.141592656589793) real a(Lorder,3,2),c(Lorder,3,4,2),r real*4 inlat,inlong,outlat,outlong,outr integer mgflag,err,order integer l,m,q real ylm,ylmval,theta,phi,x,y,z,rtheta,rphi real height real oldheigt(2) real hlevels(4) c initialize oldheigt to impossible value data oldheigt /-1.,-1./ data hlevels /0.,150.,300.,450./ c fill array c with spherical harmonic coeff.s data (c(i,1,1,1),i=1,25)/ 0.019280, 0.010959, 0.001843,-0.989591, + 0.001166, 0.012110, 0.039810, 0.015467,-0.000606,-0.000256, + -0.000982, 0.005008,-0.020183, 0.008789,-0.000889,-0.000430, + -0.000001,-0.000200, 0.000984,-0.001482,-0.009531, 0.000893, + -0.000027, 0.000124,-0.000016/ data (c(i,2,1,1),i=1,25)/ 0.013711, 0.995639, 0.004451,-0.004563, + 0.001131, 0.036807, 0.035061,-0.015005, 0.000823, 0.000304, + -0.004039, 0.005728,-0.033006,-0.005979,-0.002779,-0.000470, + 0.000063, 0.000360, 0.000776,-0.006990, 0.011303, 0.002073, + 0.000913, 0.000135, 0.000029/ data (c(i,3,1,1),i=1,25)/-0.035291,-0.073536, 1.000473, 0.043476, + 0.016315,-0.011485, 0.053083, 0.021130, 0.025971, 0.001139, + 0.002976, 0.017808,-0.001684,-0.011018, 0.002139,-0.002998, + -0.000170, 0.000100,-0.000423,-0.000542,-0.014966,-0.002059, + -0.000388,-0.000127, 0.000179/ data (c(i,1,2,1),i=1,25)/ 0.018357, 0.009213, 0.002170,-0.976504, + 0.001356, 0.011784, 0.038604, 0.014751, 0.000147,-0.000183, + -0.000883, 0.005184,-0.019307, 0.007269,-0.000931,-0.000349, + -0.000001,-0.000186, 0.000852,-0.001387,-0.009416, 0.001016, + -0.000232, 0.000091,-0.000024/ data (c(i,2,2,1),i=1,25)/ 0.012990, 0.985056, 0.004889,-0.003792, + 0.001732, 0.035455, 0.034873,-0.014480, 0.000938, 0.000228, + -0.003922, 0.005253,-0.032597,-0.005934,-0.002829,-0.000462, + 0.000038, 0.000322, 0.000548,-0.006462, 0.008353, 0.002032, + 0.000741, 0.000125, 0.000025/ data (c(i,3,2,1),i=1,25)/-0.031551,-0.064469, 1.012776, 0.052358, + 0.014111,-0.010193, 0.046556, 0.019120, 0.023080, 0.000961, + 0.002552, 0.014863,-0.017217,-0.014001, 0.001981,-0.002679, + -0.000151, 0.000101,-0.000247,-0.001455,-0.012367,-0.001380, + -0.000168,-0.000083, 0.000159/ data (c(i,1,3,1),i=1,25)/ 0.017535, 0.008299, 0.002420,-0.965697, + 0.001358, 0.011427, 0.037155, 0.014293, 0.000100,-0.000170, + -0.000797, 0.004940,-0.018514, 0.006824,-0.000926,-0.000313, + 0.000000,-0.000177, 0.000802,-0.001328,-0.008691, 0.000878, + -0.000253, 0.000078,-0.000024/ data (c(i,2,3,1),i=1,25)/ 0.012765, 0.974725, 0.004908,-0.003370, + 0.001672, 0.034194, 0.033135,-0.013955, 0.000694, 0.000187, + -0.003744, 0.004765,-0.031321,-0.005661,-0.002728,-0.000397, + 0.000035, 0.000301, 0.000517,-0.006083, 0.007747, 0.001962, + 0.000711, 0.000118, 0.000023/ data (c(i,3,3,1),i=1,25)/-0.028916,-0.059294, 1.023200, 0.061752, + 0.012918,-0.009640, 0.042250, 0.017536, 0.021127, 0.000891, + 0.002323, 0.013591,-0.029199,-0.017211, 0.001695,-0.002457, + -0.000138, 0.000100,-0.000187,-0.001413,-0.010953,-0.001058, + -0.000064,-0.000075, 0.000139/ data (c(i,1,4,1),i=1,25)/ 0.016753, 0.007485, 0.002642,-0.955238, + 0.001358, 0.011086, 0.035801, 0.013862, 0.000050,-0.000157, + -0.000718, 0.004707,-0.017769, 0.006414,-0.000922,-0.000278, + 0.000001,-0.000168, 0.000757,-0.001273,-0.008041, 0.000753, + -0.000270, 0.000067,-0.000023/ data (c(i,2,4,1),i=1,25)/ 0.012530, 0.964630, 0.004910,-0.002970, + 0.001611, 0.032991, 0.031537,-0.013457, 0.000483, 0.000150, + -0.003577, 0.004339,-0.030119,-0.005412,-0.002631,-0.000339, + 0.000033, 0.000282, 0.000490,-0.005732, 0.007205, 0.001894, + 0.000683, 0.000111, 0.000021/ data (c(i,3,4,1),i=1,25)/-0.026590,-0.054765, 1.032962, 0.070108, + 0.011895,-0.009152, 0.038477, 0.016126, 0.019449, 0.000830, + 0.002125, 0.012488,-0.040414,-0.020067, 0.001436,-0.002268, + -0.000127, 0.000099,-0.000141,-0.001355,-0.009733,-0.000786, + 0.000015,-0.000070, 0.000124/ data (c(i,1,1,2),i=1,25)/-0.016433,-0.008379,-0.004386,-1.009670, + -0.000451,-0.013717,-0.042801,-0.015660, 0.001604, 0.000517, + 0.001136,-0.006112, 0.021150,-0.008327, 0.001573, 0.000637, + -0.000014, 0.000282,-0.001095, 0.003392, 0.008625,-0.002386 , + 0.000005,-0.000155,-0.000008/ data (c(i,2,1,2),i=1,25)/-0.013804, 1.008158,-0.003596, 0.001256, + -0.000142,-0.037756,-0.029660, 0.015332, 0.000200,-0.000165, + 0.004185,-0.007108, 0.036805, 0.007977, 0.003439, 0.000620, + -0.000075,-0.000344,-0.000959, 0.010526,-0.015644,-0.003793, + -0.001329,-0.000298,-0.000042/ data (c(i,3,1,2),i=1,25)/ 0.032795, 0.071056, 0.982617,-0.042197, + -0.014513,-0.004160,-0.051034,-0.014432,-0.024027,-0.000885, + -0.001079,-0.019596, 0.017579, 0.009182,-0.000804, 0.003267, + 0.000190,-0.000014, 0.000401, 0.007867, 0.019234,-0.001521, + 0.000254,-0.000229,-0.000214/ data (c(i,1,2,2),i=1,25)/-0.014040,-0.006148,-0.005486,-1.019543, + -0.000771,-0.013548,-0.044615,-0.014943, 0.001011, 0.000451, + 0.001067,-0.006371, 0.022508,-0.008364, 0.001537, 0.000596, + -0.000011, 0.000285,-0.000970, 0.003145, 0.009593,-0.002566, + 0.000150,-0.000126,-0.000002/ data (c(i,2,2,2),i=1,25)/-0.012480, 1.017997,-0.004354,-0.000310, + -0.000590,-0.036683,-0.030536, 0.015084,-0.000344,-0.000153, + 0.004233,-0.005734, 0.036614, 0.008074, 0.003388, 0.000617, + -0.000054,-0.000324,-0.000803, 0.009513,-0.013535,-0.003687, + -0.001104,-0.000270,-0.000038/ data (c(i,3,2,2),i=1,25)/ 0.030921, 0.067449, 0.972432,-0.053186, + -0.013150,-0.003639,-0.047223,-0.013669,-0.021655,-0.000821, + -0.001035,-0.018077, 0.029247, 0.012154,-0.000871, 0.003097, + 0.000181,-0.000036, 0.000209, 0.007680, 0.016882,-0.001781, + -0.000079,-0.000238,-0.000193/ data (c(i,1,3,2),i=1,25)/-0.012001,-0.004642,-0.006475,-1.027540, + -0.000744,-0.013248,-0.045739,-0.014481, 0.000941, 0.000420, + 0.000985,-0.006345, 0.023658,-0.009124, 0.001419, 0.000568, + -0.000011, 0.000282,-0.000936, 0.002920, 0.009705,-0.002517, + 0.000183,-0.000107, 0.000000/ data (c(i,2,3,2),i=1,25)/-0.012085, 1.026935,-0.004776,-0.001409, + -0.000409,-0.035536,-0.029287, 0.014827,-0.000327,-0.000155, + 0.004205,-0.004315, 0.035758, 0.007989, 0.003252, 0.000566, + -0.000049,-0.000315,-0.000807, 0.008605,-0.013536,-0.003575, + -0.001036,-0.000244,-0.000035/ data (c(i,3,3,2),i=1,25)/ 0.029230, 0.064873, 0.962930,-0.063255, + -0.012119,-0.003165,-0.043867,-0.012826,-0.019745,-0.000795, + -0.001002,-0.017098, 0.039770, 0.014767,-0.000837, 0.002959, + 0.000174,-0.000052, 0.000075, 0.007277, 0.014957,-0.001998, + -0.000328,-0.000230,-0.000176/ data (c(i,1,4,2),i=1,25)/-0.010060,-0.003278,-0.007417,-1.035223, + -0.000717,-0.012950,-0.046770,-0.014058, 0.000864, 0.000391, + 0.000907,-0.006304, 0.024711,-0.009936, 0.001309, 0.000540, + -0.000010, 0.000279,-0.000906, 0.002712, 0.009801,-0.002461, + 0.000216,-0.000089, 0.000001/ data (c(i,2,4,2),i=1,25)/-0.011770, 1.035414,-0.005173,-0.002405, + -0.000219,-0.034423,-0.027984, 0.014579,-0.000289,-0.000155, + 0.004178,-0.002882, 0.034954, 0.007896, 0.003126, 0.000516, + -0.000045,-0.000308,-0.000812, 0.007776,-0.013613,-0.003466, + -0.000977,-0.000221,-0.000032/ data (c(i,3,4,2),i=1,25)/ 0.027660, 0.062635, 0.953620,-0.072585, + -0.011199,-0.002781,-0.040773,-0.012008,-0.018038,-0.000774, + -0.000970,-0.016244, 0.050048, 0.017101,-0.000808, 0.002842, + 0.000167,-0.000065,-0.000039, 0.006930, 0.013211,-0.002205, + -0.000539,-0.000221,-0.000162/ if (height.lt.0. .or. height.gt.2000.) then err = -2 else if (mgflag.lt.1 .or. mgflag.gt.2) then err = -4 else if ((abs(inlat).gt.90.) .or. (abs(inlong).gt.360.)) then err = -8 else if ((order .lt. 1) .or. (order.gt.4)) then err = -16 else err = 0 end if c quit if an error condition exists if (err.ne.0) return if (height.ne.oldheigt(mgflag)) then call cinterp(c(1,1,1,mgflag),hlevels,height,a(1,1,mgflag)) oldheigt(mgflag) = height end if x = 0 y = 0 z = 0 if (mgflag .eq. 2) then theta = (90.0 - inlat)*pi/180.00 phi = inlong*pi/180.00 do 200 l=0,order do 200 m=-l,l q = l*(l+1) + m + 1 ylmval = ylm(l,m,theta,phi) x = a(q,1,mgflag)*ylmval + x y = a(q,2,mgflag)*ylmval + y z = a(q,3,mgflag)*ylmval + z 200 continue r = sqrt(x*x + y*y + z*z) if (r.gt.1.1 .or. r.lt.0.9) then err = -32 end if z = z/r x = x/r y = y/r outr = r if (z.gt.1.00) then rtheta = 0.00 else if (z.lt.-1.00) then rtheta = pi else rtheta = acos(z) end if if (abs(x).lt.1.0e-8 .and. abs(y).lt.1.0e-8) then rphi = 0.0 else rphi = atan2(y,x) end if call rotatepl(rtheta,rphi,theta,phi,mgflag) outlat = 90.00 - 180.00*theta/pi outlong = 180.00*phi/pi else theta = (90.0 - inlat)*pi/180. phi = inlong*pi/180. call rotatepl(theta,phi,rtheta,rphi,mgflag) do 300 l=0,order do 300 m=-l,l q = l*(l+1) + m + 1 ylmval = ylm(l,m,rtheta,rphi) x = a(q,1,mgflag)*ylmval + x y = a(q,2,mgflag)*ylmval + y z = a(q,3,mgflag)*ylmval + z 300 continue r = sqrt(x*x + y*y + z*z) if (r.gt.1.1 .or. r.lt.0.9) then err = -32 end if z = z/r x = x/r y = y/r outr = r if (z.gt.1.00) then theta = 0.0 else if (z.lt.-1.00) then theta = pi else theta = acos(z) end if outlat = 90.00 - 180.00*theta/pi if (abs(x).lt.1.0e-8 .and. abs(y).lt.1.0e-8) then phi = 0.0 else phi = atan2(y,x) end if outlong = 180.00*phi/pi end if return end c ************************************************************************** subroutine cinterp(c,hlev,h,a) c c This subroutine interpolates the Ylm coefficients for the specified c height h from the values determined at the heights specified by c array hlev c real c(25,3,4),a(25,3) real hlev(4),h real xa(4),ya(4),x,y,dy x = h do 100 k = 1,4 xa(k) = hlev(k) 100 continue do 300 i=1,25 do 300 j=1,3 do 200 k=1,4 ya(k) = c(i,j,k) 200 continue call polint(xa,ya,4,x,y,dy) a(i,j) = y 300 continue return end c ************************************************************************** subroutine polint(xa,ya,n,x,y,dy) c c This subroutine performs a polynomial interpolation with error c estimate using Neville's algorithm. The program is taken from c the book, "numerical recipies" by Press et al., pp 80-83. c c N is the number of data points in the arrays xa and ya. The c polynomial is therefore of oreder n-1. c parameter ( nmax=10 ) real xa(n),ya(n),c(nmax),d(nmax),x,y,dy real dif,dift,ho,hp,w,den ns=1 dif = abs(x-xa(1)) c search for the element in xa closest to x do 100 i=1,n dift = abs(x-xa(i)) if (dift .lt. dif) then ns = i dif = dift end if c initialize differences c(i) = ya(i) d(i) = ya(i) 100 continue y = ya(ns) ns= ns - 1 do 400 m=1,n-1 do 200 i=1,n-m ho = xa(i) - x hp = xa(i+m) - x w = c(i+1) - d(i) den = ho - hp if (den .eq. 0.) stop 'POLINT: invalid input arrays' den = w/den d(i) = hp*den c(i) = ho*den 200 continue if (2*ns .lt. n-m) then dy = c(ns+1) else dy=d(ns) ns = ns - 1 end if y = y + dy 400 continue return end c **************************************************************************** subroutine rotatepl(intheta,inphi,outtheta,outphi,mgflag) c c This subroutine rotates geographic theta and phi into c a new coordinate system with the z axis along the magnetic c pole OR the reverse rotation. c parameter ( pi=3.141592653589793 ) real intheta,inphi,outtheta,outphi,x,y,z real thetapl,phipole,xin(3),xout(3),rmag,R(3,3) integer mgflag equivalence (x,xout(1)),(y,xout(2)),(z,xout(3)) logical first data thetapl,phipole / 0.195476876,-1.234994979 / data first /.true./ if (first) then r(1,1) = cos(thetapl)*cos(phipole) r(1,2) = cos(thetapl)*sin(phipole) r(1,3) = -sin(thetapl) r(2,1) = -sin(phipole) r(2,2) = cos(phipole) r(2,3) = 0. r(3,1) = sin(thetapl)*cos(phipole) r(3,2) = sin(thetapl)*sin(phipole) r(3,3) = cos(thetapl) first = .false. end if xin(1) = sin(intheta)*cos(inphi) xin(2) = sin(intheta)*sin(inphi) xin(3) = cos(intheta) if (mgflag .eq. 1) then do 100 j=1,3 xout(j) = 0. do 100 k = 1,3 xout(j) = xout(j) + R(j,k)*xin(k) 100 continue else if (mgflag .eq. 2) then do 200 j = 1,3 xout(j) = 0. do 200 k = 1,3 xout(j) = xout(j) + R(k,j)*xin(k) 200 continue else do 300 j=1,3 xout(j) = 0. 300 continue end if rmag = sqrt(x*x + y*y + z*z) x=x/rmag y=y/rmag z=z/rmag if (z.gt.1.0) then outtheta = 0.00 else if (z.lt.-1.0)then outtheta = pi else outtheta = acos(z) end if if (abs(x).lt.1.0e-8 .and. abs(y).lt.1.0e-8) then outphi = 0.00 else outphi = atan2(y,x) end if return end c *************************************************************************** real function ylm(l,m,theta,phi) c c This function returns the appropriate spherical harmonic c for a given latitude and longitude. This is actually c returning only the real part of Y(l,m). c integer l,m real theta,phi,plgndr if (m .lt. 0) then ylm = plgndr(l,-m,cos(theta))*sin(m*phi) else ylm = plgndr(l,m,cos(theta))*cos(m*phi) end if return end c **************************************************************************** real function plgndr(l,m,x) c c This routine computes the associated legendre polynomial P(l,m)(x) c real x,somx2,pmm,fact,pmmp1,pll integer l,m integer i,ll if (m.lt.0 .or. m.gt.l .or. abs(x).gt.1) stop 'illegal arguments' pmm = 1 if (m.gt.0) then somx2 = sqrt((1.00 - x)*(1.00 + x)) fact = 1.00 do 100 i=1,m pmm = -pmm*fact*somx2 fact = fact+2.0 100 continue end if if (l.eq.m) then plgndr=pmm else pmmp1 = x*(2*m+1)*pmm if (l.eq.m+1) then plgndr=pmmp1 else do 200 ll=m+2,l pll = (x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll 200 continue plgndr = pll end if end if return end