rotate.f
subroutine rotate (glt, gln, phir, thetar, psir, rlt, rln) c c double precision translator for drotate c c real glt, gln, phir, thetar, psir, rlt, rln real8 dglt, dgln, dphir, dthetar, dpsir, drlt, drln c dglt = glt dgln = gln dphir = phir dthetar = thetar dpsir = psir drlt = rlt drln = rln call drotate (dglt, dgln, dphir, dthetar, dpsir, drlt, drln) rlt = drlt rln = drln c return end subroutine drotate (glt, gln, phir, thetar, psir, rlt, rln) c c======================================================================= c subroutine drotate takes a geographic latitude and longitude and c finds the the equivalent latitude and longitude on a rotated grid. c when going from a geographic grid to a rotated grid, all of the c defined rotation angles given to drotate by the calling program c are positive, but when going from a rotated grid back to the c geographic grid, the calling program must reverse the angle order c (phir and psir are switched) and all of the angles made negative. c c the first rotation angle phir is defined as a rotation about the c original z axis. the second rotation angle thetar is defined as a c rotation about the new x axis. the final rotation angle psir is c defined as a rotation about the new z axis. these rotation angles c are just the Euler angles as defined in "classical mechanics" c Goldstein (1951). c c author: M. Eby e-mail [email protected] c======================================================================= c c real8 glt, gln, phir, thetar, psir, phis, thetas, rlt, rln &, gy, gx, gz, rx, ry, rz, d0, d1, d90, d180, d360, rad c c g... = geographic value c r... = rotated value c ...lt = latitude (or equivalent spherical coordinate) c ...ln = longitude (or equivalent spherical coordinate) c ...x = x coordinate c ...y = y coordinate c ...z = z coordinate c psir, thetar, phir = Euler angles defining rotation c d0 = 0.D+00 d1 = 1.D+00 d90 = 90.D+00 d180 = 180.D+00 d360 = 360.D+00 c c define rad for conversion to radians. rad = acos(-d1)/d180 c c convert latitude and longitude to spherical coordinates thetas = gln if (thetas .gt. d180) thetas = thetas - d360 if (thetas .lt. -d180) thetas = thetas + d360 thetas = thetasrad phis = (d90 - glt)rad c c translate point into Cartesian coordinates for rotation. gx = sin(phis)cos(thetas) gy = sin(phis)sin(thetas) gz = cos(phis) c c rotate the point (gx, gy, gz) about the z axis by phir then the x c axis by thetar and finally about the z axis by psir. c rx = gx(cos(psir)cos(phir) - cos(thetar)sin(phir)sin(psir)) + & gy(cos(psir)sin(phir) + cos(thetar)cos(phir)sin(psir)) + & gzsin(psir)sin(thetar) c ry = gx(-sin(psir)cos(phir) - cos(thetar)sin(phir)cos(psir)) + & gy(-sin(psir)sin(phir) + cos(thetar)cos(phir)cos(psir)) + & gz(cos(psir)sin(thetar)) c rz = gx(sin(thetar)sin(phir)) + gy(-sin(thetar)cos(phir)) + & gz(cos(thetar)) c c convert rotated point back to spherical coordinates c c check for rounding error (arccos(x): abs(x) must be .le. 1) rz = min(rz, d1) rz = max(rz, -d1) rlt = acos(rz) c if point is at a pole set rotated longitude equal to initial. if (rlt .le. d0 .or. rlt .ge. d180rad) then rln = thetas else c if rln lies between -135 and -45 or between 45 and 135 degrees c it is more accurate to use an arccos calculation. if (abs(rx/sin(rlt)) .lt. cos(45.D+00rad)) then rln = acos(max(min(rx/sin(rlt), d1), -d1)) c arccos will give rln between 0 and 180 degrees. if the point c is negative in y, rln must be equal to negative rln. if (ry .lt. d0) rln = -rln else c if rln lies between -45 and 45 or between 135 and -135 degrees c it is more accurate to use an arcsin calculation. rln = asin(max(min(ry/sin(rlt), d1), -d1)) c arcsin will give rln between -90 and 90 degrees. if the point c is negative in x, rln must be equal to 180 degrees minus rln. if (rx .lt. d0) rln = d180rad - rln endif endif c c convert back to degrees of latitude and longitude. rlt = d90 - rlt/rad rln = rln/rad if (rln .gt. d180) rln = rln - d360 if (rln .le. -d180) rln = rln + d360 return end
subroutine rotate (glt, gln, phir, thetar, psir, rlt, rln)
c
c double precision translator for drotate
c
c
real glt, gln, phir, thetar, psir, rlt, rln
real*8 dglt, dgln, dphir, dthetar, dpsir, drlt, drln
c
dglt = glt
dgln = gln
dphir = phir
dthetar = thetar
dpsir = psir
drlt = rlt
drln = rln
call drotate (dglt, dgln, dphir, dthetar, dpsir, drlt, drln)
rlt = drlt
rln = drln
c
return
end
subroutine drotate (glt, gln, phir, thetar, psir, rlt, rln)
c
c=======================================================================
c subroutine drotate takes a geographic latitude and longitude and
c finds the the equivalent latitude and longitude on a rotated grid.
c when going from a geographic grid to a rotated grid, all of the
c defined rotation angles given to drotate by the calling program
c are positive, but when going from a rotated grid back to the
c geographic grid, the calling program must reverse the angle order
c (phir and psir are switched) and all of the angles made negative.
c
c the first rotation angle phir is defined as a rotation about the
c original z axis. the second rotation angle thetar is defined as a
c rotation about the new x axis. the final rotation angle psir is
c defined as a rotation about the new z axis. these rotation angles
c are just the Euler angles as defined in "classical mechanics"
c Goldstein (1951).
c
c author: M. Eby e-mail [email protected]
c=======================================================================
c
c
real*8 glt, gln, phir, thetar, psir, phis, thetas, rlt, rln
&, gy, gx, gz, rx, ry, rz, d0, d1, d90, d180, d360, rad
c
c g... = geographic value
c r... = rotated value
c ...lt = latitude (or equivalent spherical coordinate)
c ...ln = longitude (or equivalent spherical coordinate)
c ...x = x coordinate
c ...y = y coordinate
c ...z = z coordinate
c psir, thetar, phir = Euler angles defining rotation
c
d0 = 0.D+00
d1 = 1.D+00
d90 = 90.D+00
d180 = 180.D+00
d360 = 360.D+00
c
c define rad for conversion to radians.
rad = acos(-d1)/d180
c
c convert latitude and longitude to spherical coordinates
thetas = gln
if (thetas .gt. d180) thetas = thetas - d360
if (thetas .lt. -d180) thetas = thetas + d360
thetas = thetas*rad
phis = (d90 - glt)*rad
c
c translate point into Cartesian coordinates for rotation.
gx = sin(phis)*cos(thetas)
gy = sin(phis)*sin(thetas)
gz = cos(phis)
c
c rotate the point (gx, gy, gz) about the z axis by phir then the x
c axis by thetar and finally about the z axis by psir.
c
rx = gx*(cos(psir)*cos(phir) - cos(thetar)*sin(phir)*sin(psir)) +
& gy*(cos(psir)*sin(phir) + cos(thetar)*cos(phir)*sin(psir)) +
& gz*sin(psir)*sin(thetar)
c
ry = gx*(-sin(psir)*cos(phir) - cos(thetar)*sin(phir)*cos(psir)) +
& gy*(-sin(psir)*sin(phir) + cos(thetar)*cos(phir)*cos(psir)) +
& gz*(cos(psir)*sin(thetar))
c
rz = gx*(sin(thetar)*sin(phir)) + gy*(-sin(thetar)*cos(phir)) +
& gz*(cos(thetar))
c
c convert rotated point back to spherical coordinates
c
c check for rounding error (arccos(x): abs(x) must be .le. 1)
rz = min(rz, d1)
rz = max(rz, -d1)
rlt = acos(rz)
c if point is at a pole set rotated longitude equal to initial.
if (rlt .le. d0 .or. rlt .ge. d180*rad) then
rln = thetas
else
c if rln lies between -135 and -45 or between 45 and 135 degrees
c it is more accurate to use an arccos calculation.
if (abs(rx/sin(rlt)) .lt. cos(45.D+00*rad)) then
rln = acos(max(min(rx/sin(rlt), d1), -d1))
c arccos will give rln between 0 and 180 degrees. if the point
c is negative in y, rln must be equal to negative rln.
if (ry .lt. d0) rln = -rln
else
c if rln lies between -45 and 45 or between 135 and -135 degrees
c it is more accurate to use an arcsin calculation.
rln = asin(max(min(ry/sin(rlt), d1), -d1))
c arcsin will give rln between -90 and 90 degrees. if the point
c is negative in x, rln must be equal to 180 degrees minus rln.
if (rx .lt. d0) rln = d180*rad - rln
endif
endif
c
c convert back to degrees of latitude and longitude.
rlt = d90 - rlt/rad
rln = rln/rad
if (rln .gt. d180) rln = rln - d360
if (rln .le. -d180) rln = rln + d360
return
end