anchor.diff:  
         IMPLICIT none
         REAL*8 pi,d2r
  -      PARAMETER (pi=4.0d0*atan(1.0d0),d2r=pi/1.8d2)
  +      PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
         REAL*8 lon,lat,x,y,z
   
  -      x = cos(lon*d2r)*cos(lat*d2r)
  -      y = sin(lon*d2r)*cos(lat*d2r)
  -      z = sin(lat*d2r)
  +      x = dcos(lon*d2r)*dcos(lat*d2r)
  +      y = dsin(lon*d2r)*dcos(lat*d2r)
  +      z = dsin(lat*d2r)
   
         RETURN
         END
   
   C----------------------------------------------------------------------C
   
  -      SUBROUTINE gcmdln(ifile,ofile)
  +      SUBROUTINE gcmdln(ifile,ofile,user)
   
         IMPLICIT NONE
         CHARACTER*30 ifile,ofile,tag
  -      INTEGER i,narg
  +      INTEGER i,narg,user
   
  -      ifile = 'none'
  +      user = 0
  +      ifile = 'sphfinrot.in'
         ofile = 'sphfinrot.out'
   
         narg = iargc()
  +      if (narg.eq.0) call usage()
  +
         i = 0
    9998 i = i + 1
         if (i.gt.narg) goto 9999
  
positive.code:       IMPLICIT none
  
        write(*,*)
       1 'Usage: sphfinrot -f [IFILE] -o [OFILE] -h/-?'
        write(*,*)
       1 '  -f [IFILE] (Default none) name of input file'
        write(*,*)
       1 '                            if no file name given, fault ',
       2                            'parameters prompted for manual input'
        write(*,*)
       1 '  -o [OFILE] (Default sphfinrot.out) name of output file'
        write (*,*)
       1 '  -h/-?        help'
        write (*,*) ''
        write (*,*)
       1 '  sphfinrot makes a counter-clockwise finite rotation of a',
       2     ' point around a pole'
        write (*,*) ''
        write (*,*)
       1 '    input file'
        write (*,*)
       1 '        lon_init lat_init pole_lon pole_lat angle'
        write (*,*)
       1 '         :  :'
        write (*,*)
       1 '    output file'
        write (*,*)
       1 '        lon_final lat_final'
        write (*,*)
  
negative.code:       REAL*8 olon,olat,ilon,ilat,plon,plat,angle
        REAL*8 ix,iy,iz,px,py,pz,ox,oy,oz
        REAL*8 pi,d2r,r2d
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2,r2d=1.8d2/pi)
  
        call lonlat2xyz(ix,iy,iz,ilon,ilat)
        call lonlat2xyz(px,py,pz,plon,plat)
  
        ox = px*(px*ix+py*iy+pz*iz)*(1.0d0-dcos(angle*d2r))
       1         + ix*dcos(angle*d2r)
       2         + (-pz*iy+py*iz)*dsin(angle*d2r)
        oy = py*(px*ix+py*iy+pz*iz)*(1.0d0-dcos(angle*d2r))
       1         + iy*dcos(angle*d2r)
       2         + ( pz*ix-px*iz)*dsin(angle*d2r)
        oz = pz*(px*ix+py*iy+pz*iz)*(1.0d0-dcos(angle*d2r))
       1         + iz*dcos(angle*d2r)
       2         + (-py*ix+px*iy)*dsin(angle*d2r)
  
        olon = datan2(oy,ox)*r2d
        olat = dasin(oz)*r2d
  
        RETURN
        END
  