ppf_lib_f.f

read envisat and analyis sar data from NASA

C********************************************************************
C* This version number is the SCCS version number valid only for
C* internal configuration and ontrol.
C* The oficial version number can be found in the History box below.
C*
C* SCCS ID "@(#)ppf_lib_F.F	1.7		00/04/12"
C*
C* Module      : ppf_lib_F.F
C*
C* Functions   : main program (ppf_lib example) to call:
C*
C*                    - PL_car_geo
C*                    - PL_srar_cs
C*                    - PL_geo_distance
C*                    - tmjd
C*                    - emjd 
C*                    - pmjd 
C*                    - tadd 
C*                    - tsub 
C*                    - PL_sun
C*                    - PL_moon
C*
C* Purpose     : Example program to show the way to call the 
C*               previous functions
C*
C*               +-------------------------------------------------+
C* History     : |Version |  Date  |   Name    | Change            |
C*               |-------------------------------------------------|
C*               |   4.1  |15/01/97| GMV, S.A. | First release     |
C*               |   4.2  |30/04/97| GMV, S.A. | Second release    |
C*               |   4.3  |18/05/98| GMV, S.A. | Third release     |
C*               |   4.4  |13/10/98| GMV, S.A. | Fourth release    |
C*               |   4.5  |25/05/99| GMV, S.A. | Fifth release     |
C*               |   4.6  |14/04/00| GMV, S.A. | Sixth release     |
C*               |   4.7  |22/06/01| GMV, S.A. | Seventh release   |
C*               |   4.8  |31/07/01| GMV, S.A. | Eigth release     |
C*               |   4.9  |22/10/01| DEIMOS    | Ninth release     |
C*               |   5.0  |18/01/02| DEIMOS    | Tenth release     |
C*               |   5.1  |25/11/02| DEIMOS    | Eleventh release  |
C*               |   5.2  |26/05/03| DEIMOS    | Twelfth release   |
C*               |   5.3  |13/12/04| DEIMOS    | Thirteenth release|
C*               |  5.3.1 |15/02/05|ESA/EOP-PES| Mac OS X / Linux  |
C*               |        |        |           | / IRIX Porting    |
C*               |   5.4  |17/05/05| DEIMOS    | Fifteenth release |
C*               +-------------------------------------------------+
C*
C******************************************************************/

      PROGRAM MAIN

      IMPLICIT NONE

#include "ppf_lib.inc"

      INTEGER*4 STATUS

      INTEGER*4 MJDT(4), MJDT_2(4), MJDT_3(4)
      REAL*8 MJDP(2), MJDP_2(2), MJDP_3(2)

      CHARACTER UTCE*28, DUT1E*9, UTCE_2*28
      CHARACTER DUT1E_2*9, UTCE_3*28, DUT1E_3*9

      INTEGER*4 MJDT_A1(4), MJDT_ADD(4), MJDT_A3(4), MJDT_SUB(4)

      REAL*8 UT1(2), RSUN(3), RDSUN(3)

      REAL*8 UT1_2(2), RMOON(3), RDMOON(3)

      INTEGER*4 PLANET_ID  

      REAL*8 R_PLANET(3), RD_PLANET(3)

      REAL*8 RA0, DEC0, MU_RA, MU_DEC, RAD_VEL, PAR, RA, DEC

      INTEGER*4 CS_IN_T, CS_OUT_T
      REAL*8 TIME_T(2), R_T(3), RD_T(3), R2D_T(3)
      REAL*8 R_OUT_T(3), RD_OUT_T(3), R2D_OUT_T(3)

      REAL*8 RR_T(3), RRD_T(3), RR2D_T(3)
      REAL*8 LON_T, LAT_T, H_T, LOND_T, LATD_T, HD_T

      REAL*8 AOCS_T(3), MISP_T(3), MISPD_T(3)
      REAL*8 XS_T(3), YS_T(3), ZS_T(3)
      REAL*8 XSD_T(3), YSD_T(3), ZSD_T(3)
      REAL*8 ANGOUT_T(3), ANGOUTD_T(3)
      REAL*8 LON1_T, LAT1_T, LON2_T, LAT2_T
      REAL*8 D_T, AZ1_T, AZ2_T
   
      INTEGER*4 EXT_STATUS

C**** EXT_STATUS is the (extended) status flag ****
 
      INTEGER*4 FUNC_ID
      INTEGER*4	N
      CHARACTER*256 MSG(64)

      REAL*8 DUMMY_D
      INTEGER*4 DUMMY_L

      DUMMY_D = 0.0
      DUMMY_L = 0

C*****************************************************************
C*                                                               *
C     LOGGING FEATURE
C*                                                               *
C*****************************************************************
 
      N=2
      MSG(1) = "This is the first user's log message\0"
      MSG(2) = "This is the second user's log message\0"
      STATUS =  PL_PRINT_MSG(N,MSG)

C*****************************************************************
C*                                                               *
C* Change of coordinates from EARTH_FIXED to TRUE_DATE           *
C*                                                               *
C*   CS_IN_T is the initial coordinate  system                   *
C*   CS_OUT_T is the final coordinate system                     * 
C*   CS_IN_T and CS_OUT_T in accordance with the C enumeration   *
C*                                                               *
C*   TIME_T(0) is the UTC time in MJD2000 [days]                 *
C*   TIME_T(1) is the Delta UT1 [s]                              *
C*                                                               *
C*   R_T is the position vector (x,y,z) [m]                      *        
C*   RD_T is the velocity vector (x,y,z) [m/s]                   *
C*   R2D_T is the acceleration vector (x,y,z) [m/(s*s)]          *
C*   R_OUT_T is the position vector (x,y,z) [m]                  *        
C*   RD_OUT_T is the velocity vector (x,y,z) [m/s]               *
C*   R2D_OUT_T is the acceleration vector (x,y,z) [m/(s*s)]      *
C*                                                               *
C*****************************************************************
         
C  1993-04-11    00:00:  .00    
      TIME_T(1)= -2456.0
      TIME_T(2)= 0.D0
 
      CS_IN_T = 6
      CS_OUT_T = 5

C***  CS_IN_T=6 means Earth Fixed Coordinate System  ****
C***  CS_OUT_T=5 means True of Date Coordinate System ***

      R_T(1) = 4859964.138
      R_T(2) = -5265612.059
      R_T(3) = 0000.002
      RD_T(1) = -1203.303801
      RD_T(2) = -1098.845511
      RD_T(3) = 7377.224410
      R2D_T(1) = 0.D0
      R2D_T(2) = 0.D0
      R2D_T(3) = 0.D0

      N = 1
      MSG(1) = "\n\nPL_CHANGE_SV_CS\n\0"
      STATUS =  PL_PRINT_MSG(N,MSG) 
      N = N+1     
                       
      EXT_STATUS=PL_CHANGE_SV_CS(TIME_T, CS_IN_T, CS_OUT_T, R_T, RD_T,
     &                     R2D_T, R_OUT_T, RD_OUT_T, R2D_OUT_T)
      
      IF (EXT_STATUS .NE. 0) THEN           
          FUNC_ID = 0
          STATUS = PL_VECTOR_MSG(FUNC_ID, EXT_STATUS, N, MSG)
          STATUS = PL_PRINT_MSG(N, MSG)

          IF (EXT_STATUS.LE. -1) THEN
          STOP 'PL_CHANGE_SV_CS failed'
          ENDIF

      ENDIF 

      WRITE(MSG(1),400) 'R_OUT_T(1)', R_OUT_T(1), '\0' 
      WRITE(MSG(2),400) 'R_OUT_T(2)', R_OUT_T(2), '\0' 
      WRITE(MSG(3),400) 'R_OUT_T(3)', R_OUT_T(3), '\0'
      WRITE(MSG(4),400) 'RD_OUT_T(1)', RD_OUT_T(1), '\0'
      WRITE(MSG(5),400) 'RD_OUT_T(2)', RD_OUT_T(2), '\0'
      WRITE(MSG(6),400) 'RD_OUT_T(3)', RD_OUT_T(3), '\0'
      WRITE(MSG(7),400) 'R2D_OUT_T(1)', R2D_OUT_T(1), '\0'
      WRITE(MSG(8),400) 'R2D_OUT_T(2)', R2D_OUT_T(2), '\0'
      WRITE(MSG(9),400) 'R2D_OUT_T(3)', R2D_OUT_T(3), '\0'
      STATUS =  PL_PRINT_MSG(9,MSG)
    
C*******************************************************/
C*                                                     */
C* Change between cartesian and geodetic coordinates   */
C*                                                     */
C*     RR_T is the position vector in Earth            */
C*     fixed coordinate system (x,y,z) [m]             */   
C*     RRD_T is the velocity vector in Earth           */
C*     fixed coordinate system (Vx,Vy,Vz) [m]          */ 
C*     RR2D_T is the acceleration vector in Earth      */
C*     fixed coordinate system (Vx,Vy,Vz) [m]          */ 
C*                                                     */
C*     LON_LAT is the geocentric longitude [deg]       */
C*     LAT_T is the geodetic latitude [deg]            */
C*     H_T is the geodetic altitude [m]                */ 
C*     LOND_T is the geocentric longitude rate [deg/s] */
C*     LATD_T is the geodetic latitude rate [deg/s]    */        
C*     HD_T is the geodetic altitude rate [m/s]        */
C*                                                     */ 
C*******************************************************/

      RR_T(1) = 4859964.138
      RR_T(2) = -5265612.059
      RR_T(3) = 0000.002
      RRD_T(1) = -1203.303801
      RRD_T(2) = -1098.845511
      RRD_T(3) = 7377.224410

C Change from cartesian to geodetic coordinates 
C************************************************

      N=1
      MSG(1) = "\n\nPL_CAR_GEO\n\0"
      STATUS =  PL_PRINT_MSG(N,MSG)
      N = N+1 

      EXT_STATUS=PL_CAR_GEO( RR_T, RRD_T, LON_T, LAT_T, H_T, LOND_T,
     &                  LATD_T, HD_T)

      IF (EXT_STATUS .NE. 0) THEN
          FUNC_ID = 1
          STATUS = PL_VECTOR_MSG(FUNC_ID, EXT_STATUS, N, MSG)
          STATUS = PL_PRINT_MSG(N, MSG)

          IF (EXT_STATUS.LE. -1) THEN
          STOP 'PL_CAR_GEO failed'
          ENDIF
      ENDIF 

      WRITE(MSG(1),400) 'LON_T', LON_T, '\0' 
      WRITE(MSG(2),400) 'LAT_T', LAT_T, '\0' 
      WRITE(MSG(3),400) 'H_T', H_T , '\0'
      WRITE(MSG(4),400) 'LOND_T', LOND_T , '\0'
      WRITE(MSG(5),400) 'LATD_T', LATD_T , '\0'
      WRITE(MSG(6),400) 'HD_T', HD_T , '\0'
      STATUS =  PL_PRINT_MSG(6,MSG)
      
C Change from geodetic to cartesian coordinates 
C************************************************

      LON_T =  4.17
      LAT_T =  45.0
      H_T = 400.0
      LOND_T = 0.D0
      LATD_T = 0.D0
      HD_T =  0.D0

      N=1
      MSG(N) = "\n\nPL_GEO_CAR\n\0"
      STATUS =  PL_PRINT_MSG(N,MSG)

      EXT_STATUS=PL_GEO_CAR( LON_T, LAT_T, H_T, LOND_T,
     &                  LATD_T, HD_T, RR_T, RRD_T)

      IF (EXT_STATUS .NE. 0) THEN
          FUNC_ID = 2
          STATUS = PL_VECTOR_MSG(FUNC_ID, EXT_STATUS, N, MSG)
          STATUS = PL_PRINT_MSG(N, MSG)

          IF (EXT_STATUS.LE. -1) THEN
          STOP 'PL_GEO_CAR failed'
          ENDIF
      ENDIF 

      WRITE(MSG(1),400) 'RR_T(1)', RR_T(1), '\0'
      WRITE(MSG(2),400) 'RR_T(2)', RR_T(2), '\0'
      WRITE(MSG(3),400) 'RR_T(3)', RR_T(3), '\0'
      WRITE(MSG(4),400) 'RRD_T(1)', RRD_T(1), '\0'
      WRITE(MSG(5),400) 'RRD_T(2)', RRD_T(2), '\0'
      WRITE(MSG(6),400) 'RRD_T(3)', RRD_T(3), '\0'

      STATUS =  PL_PRINT_MSG(6,MSG)

      
C********************************************************
C*                                                      *
C* Calculate the Satellite Relative Actual Reference    * 
C*                                                      * 
C*     AOCS_T(1) is the Cx AOCS parameter (roll) [deg]  *
C*     AOCS_T(2) is the Cy AOCS parameter (pitch) [deg] *
C*     AOCS_T(3) is the Cz AOCS parameter (yaw) [deg]   *
C*                                                      *
C*     MISP_T is the vector of the mispointing angles   *
C*            between SRAR and SRR coordinate systems   *
C*            (x,y,z) [deg]                             *
C*     MISPD_T is the vector of the mispointing rates   *
C*            between SRAR and SRR coordinate systems   *
C*            (x,y,z) [deg/s]                           *
C*                                                      *
C*     XS_T is a vector containing the cartesian        *
C*            coordinates of the unitary direction      *
C*            vector in the x-axis of the SRAR          *
C*            coordinate system expressed in the Earth  *
C*            fixed coordinate system (x,y,z) (without  *
C*            units)                                    *
C*     YS_T is a vector containing the cartesian        *
C*            coordinates of the unitary direction      *
C*            vector in the y-axis of the SRAR          *
C*            coordinate system expressed in the Earth  *
C*            fixed coordinate system (x,y,z) (without  *
C*            units)                                    *
C*     ZS_T is a vector containing the cartesian        *
C*            coordinates of the unitary direction      *
C*            vector in the z-axis of the SRAR          *
C*            coordinate system expressed in the Earth  *
C*            fixed coordinate system (x,y,z) (without  *
C*            units)                                    *
C*     XSD_T is a vector containing the cartesian       *
C*            coordinates of the first time derivative  *
C*            of the unitary direction vector of the    *
C*            x-axis in the SRAR coordinate system      *
C*            expressed in the Earth fixed coordinate   *
C*            system (x,y,z) [1/s]                      *
C*     YSD_T is a vector containing the cartesian       *
C*            coordinates of the first time derivative  *
C*            of the unitary direction vector of the    *
C*            y-axis in the SRAR coordinate system      *
C*            expressed in the Earth fixed coordinate   *
C*            system (x,y,z) [1/s]                      *
C*     ZSD_T is a vector containing the cartesian       *
C*            coordinates of the first time derivative  *
C*            of the unitary direction vector of the    *
C*            z-axis in the SRAR coordinate system      *
C*            expressed in the Earth fixed coordinate   *
C*            system (x,y,z) [1/s]                      *
C*                                                      *
C*     ANGOUT_T is a vector containing the mispointing  *
C*            angles between the SRAR and the SR        *
C*            coordinate systems (roll, pitch, yaw)     *
C*            [deg]                                     *  
C*     ANGOUTD_T is a vector containing the mispointing *
C*            rates between the SRAR and the SR         *
C*            coordinate systems (roll, pitch, yaw)     *
C*            [deg/s]                                   *      
C*                                                      *  
C********************************************************
                                             
      RR_T(1) =  4859964.138
      RR_T(2) =  -5265612.059
      RR_T(3) =  0000.002

      RRD_T(1) =  -1203.303801
      RRD_T(2) =  -1098.845511
      RRD_T(3) =  7377.224410                                      

      RR2D_T(1) =  0.D0
      RR2D_T(2) =  0.D0
      RR2D_T(3) =  0.D0

      AOCS_T(1) =  -0.1672
      AOCS_T(2) =  0.0501
      AOCS_T(3) =  3.9284 

      MISP_T(1) =  0.D0
      MISP_T(2) =  0.D0
      MISP_T(3) =  0.D0                   
   
      MISPD_T(1) =  0.D0
      MISPD_T(2) =  0.D0
      MISPD_T(3) =  0.D0

      N=1
      MSG(N) = "\n\nPL_SRAR_CS\n\0"
      STATUS =  PL_PRINT_MSG(N,MSG)
                                                         
      EXT_STATUS=PL_SRAR_CS(RR_T, RRD_T,RR2D_T, AOCS_T, MISP_T, MISPD_T,
     &        XS_T,YS_T, ZS_T, XSD_T, YSD_T, ZSD_T, ANGOUT_T, ANGOUTD_T)

      IF (EXT_STATUS .NE. 0) THEN
          FUNC_ID = 3
          STATUS = PL_VECTOR_MSG(FUNC_ID, EXT_STATUS, N, MSG)
          STATUS = PL_PRINT_MSG(N, MSG)

          IF (EXT_STATUS.LE. -1) THEN
          STOP 'PL_SRAR_CS failed'
          ENDIF
      ENDIF 

      WRITE(MSG(1),400) 'XS_T(1)', XS_T(1), '\0'
      WRITE(MSG(2),400) 'XS_T(2)', XS_T(2), '\0'
      WRITE(MSG(3),400) 'XS_T(3)', XS_T(3), '\0'
                                            
      WRITE(MSG(4),400) 'YS_T(1)', YS_T(1), '\0'
      WRITE(MSG(5),400) 'YS_T(2)', YS_T(2), '\0'
      WRITE(MSG(6),400) 'YS_T(3)', YS_T(3), '\0'
                                            
      WRITE(MSG(7),400) 'ZS_T(1)', ZS_T(1), '\0'
      WRITE(MSG(8),400) 'ZS_T(2)', ZS_T(2), '\0'
      WRITE(MSG(9),400) 'ZS_T(3)', ZS_T(3), '\0'

      WRITE(MSG(10),400) 'XSD_T(1)', XSD_T(1), '\0'
      WRITE(MSG(11),400) 'XSD_T(2)', XSD_T(2), '\0'
      WRITE(MSG(12),400) 'XSD_T(3)', XSD_T(3), '\0'
                                            
      WRITE(MSG(13),400) 'YSD_T(1)', YSD_T(1), '\0'
      WRITE(MSG(14),400) 'YSD_T(2)', YSD_T(2), '\0'
      WRITE(MSG(15),400) 'YSD_T(3)', YSD_T(3), '\0'
                                            
      WRITE(MSG(16),400) 'ZSD_T(1)', ZSD_T(1), '\0'
      WRITE(MSG(17),400) 'ZSD_T(2)', ZSD_T(2), '\0'
      WRITE(MSG(18),400) 'ZSD_T(3)', ZSD_T(3), '\0'

      WRITE(MSG(19),400) 'ANGOUT_T(1)', ANGOUT_T(1), '\0'
      WRITE(MSG(20),400) 'ANGOUT_T(2)', ANGOUT_T(2), '\0'
      WRITE(MSG(21),400) 'ANGOUT_T(3)', ANGOUT_T(3), '\0'

      WRITE(MSG(22),400) 'ANGOUTD_T(1)', ANGOUTD_T(1), '\0'
      WRITE(MSG(23),400) 'ANGOUTD_T(2)', ANGOUTD_T(2), '\0'
      WRITE(MSG(24),400) 'ANGOUTD_T(3)', ANGOUTD_T(3), '\0'