/********************************************************************************
*  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 ;