c prog to test formulas for rotation of axes
c written by leo mcnamara 30-jul-90
c----- last modified 31-jul-90
c
c
data pi,pid2,degs,rads / 3.141592654,1.570796327,57.29577951,
a 0.017453292 /
c
print*,' Prog to find new coords of given point in a rotated'
print*,' coordinate system (e.g. geogmagnetic from geogc)'
print*,' Enter coords of north pole of new coord system'
print*,' as lat & lon in old coordinate system - 2f5.0'
read(5,500) pole_lat,pole_lon
500 format(2f5.0)
10 continue
print*,' Enter coords of given point, in old coord system'
read(5,500) slat,slon
if(slat+slon.eq.0.) stop ' Normal exit'
WRITE(6,103) SLAT,SLON
103 FORMAT(/5X,'ORIGINAL LATITUDE & LONGITUDE ARE',2F6.1//)
c
print*,' Enter +1 for geog->geom; -1 for geom->geog'
read (*,*) iflag
c
if(iflag.eq.+1) geog_lat = slat
if(iflag.eq.+1) geog_lon = slon
if(iflag.eq.-1) geom_lat = slat
if(iflag.eq.-1) geom_lon = slon
call davies
a (pole_lat,pole_lon,geog_lat,geog_lon,geom_lat,geom_lon)
print*,' Davies gives ', geog_lat,geog_lon,geom_lat,geom_lon
c
if(iflag.eq.+1) geog_lat = slat
if(iflag.eq.+1) geog_lon = slon
if(iflag.eq.-1) geom_lat = slat
if(iflag.eq.-1) geom_lon = slon
call gg_gm
a (iflag,pole_lat,pole_lon,geog_lat,geog_lon,geom_lat,geom_lon)
print*,' GG_GM gives ',
a iflag,geog_lat,geog_lon,geom_lat,geom_lon
c
c this next set of coding comes from Mike Jones prog NITIAL
plat = pole_lat * rads
plon = pole_lon * rads
c
SP=SIN(PLAT)
CP=SIN(PID2-PLAT)
C
C CALCULATE COORDS IN NEW COORD SYSTEM
TLAT=SLAT*RADS
TLON=SLON*RADS
SDPH=SIN(TLON-PLON)
CDPH=SIN(PID2-(TLON-PLON))
SL=SIN(TLAT)
CL=SIN(PID2-TLAT)
THETA=Asin(CDPH*CP*CL+SP*SL)
PHI=ATAN2(SDPH*CL,CDPH*SP*CL-CP*SL)
C
C CONVERT FROM CO-LATITUDE
cccccc THETA=PID2-THETA
THETA=THETA*DEGS
C CONVERT TO DEGREES EAST LONGITUDE
PHI=PHI+2.*PI
IF(PHI.GT.2.*PI) PHI=PHI-2.*PI
PHI=PHI*DEGS
c
WRITE(6,102) THETA,PHI
102 FORMAT(//5X,'JONES LATITUDE & LONGITUDE ARE',2F6.1///)
GO TO 10
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine davies (pole_lat,pole_lon,sslat,sslon,rlat,rlon)
pi = 3.141592654
degrad = pi / 180.
radeg = 180. / pi
plat = pole_lat * degrad
plon = pole_lon * degrad
flat = sslat * degrad
flon = sslon * degrad
c eqn 1.42
sinphi = sin(flat)*sin(plat) +
a cos(flat)*cos(plat)*cos(flon-plon)
phi = asin(sinphi)
c eqn 1.43
sinlam = cos(flat)*sin(flon-plon)/cos(phi)
flam = asin(sinlam)
rlat = phi * radeg
rlon = flam * radeg
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine gg_gm
a (flag,pole_lat,pole_lon,gglat,gglon,gmlat,gmlon)
c
c s/r to convert between geographic and geomagnetic coords (both ways)
c can also be used to convert between any other coord systems
c written by leo mcnamara (from the IRI s/r GGM) 30-jul-90
c----- last modified 31-jul-90
c
real lat
integer flag
c flag = 1 for geographic to geomagnetic conversion (gg -> gm)
c flag = -1 for geomagnetic to geographic conversion
c pole_lat is the geogc latitude of the north geom pole
c pole_lon is the geogc longitude of the north geom pole
c
data pi,radeg,degrad / 3.141592654,57.29577951,0.017453292 /
c
cp = cos((90.0 - pole_lat) * degrad)
sp = sin((90.0 - pole_lat) * degrad)
c
if(flag.eq.1) slat = gglat
if(flag.eq.-1) slat = gmlat
if(flag.eq.1) slon = gglon +(360.0 - pole_lon)
if(flag.eq.-1) slon = gmlon
c
cb = cos(slat*degrad)
sb = sin(slat*degrad)
cl = cos(slon*degrad)
sl = sin(slon*degrad)
sbb = sb * cp + flag * cb * cl * sp
if(abs(sbb).gt.1.0) sbb = sign(1.0,sbb)
lat = asin(sbb)
c
if(flag.eq.1) gmlat = lat * radeg
if(flag.eq.-1) gglat = lat * radeg
c
cg = cos(lat)
cll = (-flag * sb * sp + cb * cl * cp) / cg
if(cll.gt.1.) cll = sign(1.0,cll)
flon = acos(cll)
sll = cb * sl / cg
if(sll.lt.0.) flon = 2.0 * pi - flon
c
if(flag.eq.1) gmlon = flon * radeg
if(flag.eq.-1) gglon = flon * radeg - (360.0 - pole_lon)
c
if(gmlon.gt.360.) gmlon = gmlon - 360.
if(gmlon.le.0. ) gmlon = gmlon + 360.
if(gglon.gt.360.) gglon = gglon - 360.
if(gglon.le.0. ) gglon = gglon + 360.
c
return
end