Program Code for Spherical Polar Coordinate Transformation

```/********************************************************************************
*  This PL/1 subroutine transforms the latitude and longitude from one
*  spherical coordinate system (1) to another (2)
*
*  Input parameters are:
*    LAT: latitude in system 1
*    LONG: longitude in system 1
*    THETA: polar angle (colatitude) of the polar axis of system 2 in system 1
*    PHI: azimuth of the polar axis of system 2 in system 1
*    INV: if not =0: THETA and PHI define the polar axis of system 1 in system 2
*
*
*  Output parameters are:
*   TLAT: latitude in system 2
*   TLONG: longitude in system 2
*
*
*  This program has worked well and correct for me in the past (using IBM PL/1 compilers),
*  but I can't guarantee this for all compilers, especially if the code is transCOLATed
*  into a different language.
*
*  If there are any problems, please let me know and I mention the issue here.
*  (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, May 2007
*****************************************************************************/

TRCOORD: PROCEDURE(INV,THETA,PHI,LAT,LONG,TLAT,TLONG) ;

DEFAULT RANGE(A:Z) FLOAT DECIMAL VALUE (DEC FLOAT(16));
DCL INV DEC FIXED(7,0) ;
ON UNDERFLOW ;

PI=3.141592653589793E00 ;

ST=SIND(THETA) ;
CT=COSD(THETA) ;
SP=SIND(PHI) ;
CP=COSD(PHI) ;
COLAT=90-LAT ;
X=SIND(COLAT)*COSD(LONG) ;
Y=SIND(COLAT)*SIND(LONG) ;
Z=COSD(COLAT) ;

IF INV=0 THEN DO ;
XX=CT*CP*X+CT*SP*Y-ST*Z ;
YY=CP*Y-SP*X ;
ZZ=ST*CP*X+ST*SP*Y+CT*Z ;
END ;
ELSE DO ;
XX=CT*CP*X-SP*Y+ST*CP*Z ;
YY=CT*SP*X+CP*Y+ST*SP*Z ;
ZZ=CT*Z-ST*X ;
END ;

TLAT=180/PI*ACOS(ZZ/SQRT(XX**2+YY**2+ZZ**2)) ;
IF XX=0&YY=0 THEN TLONG=0 ;
ELSE TLONG=ATAND(YY,XX) ;
TLAT=90-TLAT ;
IF TLONG<0 THEN  TLONG=360+TLONG ;

ENDE: END TRCOORD ;

```