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