* Tue May 17 13:46:05 EDT 1994 /plas1/h1/mzh/src/test phase3.f
c
c  Subroutine of MJSANL to accomodate trajectory information
c  reading from short cruise SEDR tapes or files.
c
c  The code was modified from phase2.f originally written
c  by R. Selesnick at 6/17/86.
c
c  Modifications done by: Adam Szabo   2/7/94
c
c ************************************************************
c *                                                          *
c *   From my notes of 7/23/87:                              *
c *   Spacecraft clock starts (as noted in comments in       *
c *   subroutine EDRRED):                                    *
c *                                                          *
c *   Voyager 1    1977-248/0433:29.766                      *
c *   Voyager 2    1977-232/0654:24.439                      *
c *                                                          *
c ************************************************************
c
c
      subroutine phase3(jtb,ieod)
c
c  Pathed variables:
c    jtb(6) {int*2}: Requested time array
c     jtb(1): Year in 4 digit format (i.e. 1981)
c     jtb(2): Day of year (January 1st is day 1)
c     jtb(3): Hour
c     jtb(4): Minutes
c     jtb(5): Seconds
c     jtb(6): mSeconds
c
c  The routine returns the standard read error flag ieod:
c     ieod {int*4}: 0 - OK; -1 - end of file; -2 - read error
c
c  The routine also fills the common blocks SPEC, MID, and SCID.
c  For details of the filled values of SPEC, see the comments
c  of xcru further down in this module.
c  The common block MID contains the variable IPLNT:
c    it is set to 0 in this routine to indicate cruise mode
c  The common block SCID contains the varibale VOY:
c    it will be set to 453. for Voyager 1
c                   or 372. for Voyager 2
c
c  Navigation block data variables
c
      dimension rnav(56),rnav1(56),rnav2(56),rnavt1(56),rnavt2(56)
c
c  Pointing block data variables
c
      dimension rpt(9),rpt1(9),rpt2(9),rptt1(9),rptt2(9)
c
c  Time arrays
c
      dimension itb(6),itn(6),itp(8)
c
c  Decimal hours from 1977. They need to be double precision
c  because the large values involved which are very close to each other.
c
      real*8 rhb,rhn1,rhn2,rhnt1,rhnt2,rhp1,rhp2,rhpt1,rhpt2
c
c  VGRTRJ is the single large array into which a complete short SEDR
c  reacord will be read. However, not all of its contents are reals.
c  The integer values are converted by equivalencing. For details of
c  the short SEDR record structure see Voyager Memo #180 - Revised.
c
      real*4 vgrtrj(79)
      equivalence (itn(1),vgrtrj(1)), (itp(1),vgrtrj(63))
      equivalence (rnavt2(1),vgrtrj(7)), (rptt2(1),vgrtrj(71))
      equivalence (mod16,itp(7)), (mod60,itp(8))
c
c  The requested time array JTB is passed in the calling sequence
c  as an integer*2 array.
c
      integer*2 jtb(6)
c
c  Definition of logical flags and character variables are
c  at the data section.
c

      logical ptl,navl,first,rew
      character*2 scch 
      character*2 vv1, vv2
      character*4 unit
c
c  The planet ID will be set to zero to indicate cruise mode.
c
      common/mid/iplnt 
c
c  S/C ID will be set: Voyager 1 = 453. ; Voyager 2 = 372.
c
      common/scid/voy
c
c  Initialize the temporary navigation and pointing decimal
c  time variables to zero.
c
      data rhnt2/0.d0/, rhpt2/0.d0/, rhnt1/0.d0/, rhpt1/0.d0/
c
c  Set the short SEDR file unit name as required by treadbat
c
      data unit/'ft33'/
c
c  Initialize flags signaling that pointing and navigation data
c  was already read false.
c
      data ptl/.false./, navl/.false./
c
c  Initialize flag signaling first time excution of routine
c
      data first/.true./
c
c  Initialize rewind flag false
c
      data rew/.false./
c
c  Set possible S/C names. In the MJSANL calling script, the
c  environment variable SC has to be set to either 'v1' or
c  'v2' for correct identification.
c
      data vv1 /'v1'/, vv2 /'v2'/
c
c============================================================
c      BEGIN EXECUTION OF THE ROUTINE
c============================================================
c
c  Set iplnt to 0 for cruise mission phase
c  Passed back via common block "mid"
c
      if (first) iplnt=0
c
c  Convert the requested time from int*2 to int*4
c
      do 10 i=1,6
 10      itb(i)=jtb(i)
c   
c  Convert the requested time to decimal hours from 1977 day 1
c  dechr is currently in the module phase2.f.
c
      call dechr(itb,rhb)
c
c  INPUT: itb(6) {int*4}: Usual time array.
c  OUTPUT: rhb {real*8}: Decimal hours from 1977 day 1
c
c  Check if the beginning of the short SEDR file is more than
c  10 days ahead of the requested time. If so, terminate.
c  This check is only executed after the initialization or
c  after rewinding.
c
 15   if ((.not. first) .and. (rhnt1 .eq. 0.d0)) then
         if (rhnt2 .gt. rhb+240.d0) then
            write(6,*)'ERROR! The first navigation record time in the'
            write(6,*)' SEDR file is more than 10 days ahead of the'
            write(6,*)' requested time.'
            write(6,*)
            write(6,*)'**** TERMINATING ****'
            write(6,*)
            stop
         endif
         if (rhpt2 .gt. rhb+240.d0) then
            write(6,*)'ERROR! The first pointing record time in the'
            write(6,*)' SEDR file is more than 10 days ahead of the' 
            write(6,*)' requested time.'    
            write(6,*) 
            write(6,*)'**** TERMINATING ****' 
            write(6,*) 
            stop 
         endif
      endif
c
c  Compare requested time with previous lower boundary of time
c  interval. Rewind if necessary.
c
      if (rhb.lt.rhnt1 .or. rhb.lt.rhpt1) rew=.true.
      if(rew) call tclosebat(unit)
c
c  Send through initialization only first time around or if rewinding
c
      if(.not.rew .and. .not.first) goto 49
c
c  Get spacecraft ID from environmental variable SC
c  (rather than from header record as in phase2).
c  Store S/C id in common block SCID.
c
      call getenv('SC',scch)
      if(scch .eq. vv1 ) voy=453. 
      if(scch .eq. vv2 ) voy=372. 
c
c  Write out S/C name only first time around
c
      if (first) then
         if (voy .eq. 453.) then
            write(6,*)'SEDR S/C ID: Voyager 1'
         else
            write(6,*)'SEDR S/C ID: Voyager 2'
         endif
c
c  Turn the initialization flag 'first' off
c
         first = .false.
      endif
c
c  Note if rewinding and reset time
c
      if (rew) then
         write(6,*)'Rewinding Short SEDR file'
         rhnt2 = 0.d0
         rhpt2 = 0.d0
      endif
c
c  Read records to get one on each side of rhb
c          - AND -
c  Store previous values
c
 20   rhnt1=rhnt2
      rhpt1=rhpt2
      do 30 i=1,56
 30      rnavt1(i)=rnavt2(i)
      do 40 i=1,9
 40      rptt1(i)=rptt2(i)
c
c  Read next values
c
c  Set the maxmimum number of bytes read by treadbat to length
c  of the short SEDR record. It is stored in the variable LEN.
c
      len = 4*79
      call treadbat(unit, vgrtrj, len)
c
c  TREADBAT: Routine to read from unit the next record upto
c   a maximum of len bytes into the array vgrtrj.
c   UNIT {char*4}: Name of the unit (i.e. "ft33")
c   LEN {int*4}: Maximum # of bytes allowed to be read
c   VGRTRJ(79) {real*4}: Short SEDR record will be read into this
c       array. Integer elements have to be converted before use.
c   Since treadbat is a c routine residing in the module bat.c,
c   fortran will actually call treadbat_ with a 4th variable passed.
c   This 4th variable is an integer containing the length of the
c   string UNIT. It will be received in the local variable IL.
c
      call wait(unit, len, ieod)
c
c  WAIT: Routine in the module vgrtrd.f to return the # of bytes
c   most recently read on UNIT and placed in LEN. Also the error
c   flag is returned into IEOD.
c   UNIT {char*4}: The name of the unit. The same format as for
c                  TREADBAT.
c   LEN {int*4}:   The # of bytes read most recently from UNIT.
c   IEOD {int*4}:  Error flag.   0 = Ok.
c                               -1 = End of file.
c                               -2 = End of tape.
c                               -3 = Read error.
c
c  Check for error conditions on read.
c
      if (len .le. 0) then
         write(6,*)'End of file on mission SEDR file. Terminating.'
         stop
      endif
      if ((len .ne. 316) .or. (ieod .ne. 0)) then
         write(6,*)'End of tape or read error on mission SEDR file.'
         stop
      endif
c
c  Convert to decimal hours the navigation and pointing times
c
      call dechr(itn,rhnt2)
      call dechr(itp,rhpt2)
c
c  If rewinding then go back and check if there was early
c  enough time on the SEDR file to go back to. Also turn off
c  the rewind flag.
c
      if (rew) then
         rew = .false.
         go to 15
      endif
c
c  If nav record is not in range keep on looping
c
 49   if(.not.(rhnt1.lt.rhb .and. rhb.le.rhnt2)) goto 60
c
c  If nav record is in range store data in final variables
c
         rhn1=rhnt1
         rhn2=rhnt2
         do 50 i=1,56
            rnav1(i)=rnavt1(i)
 50         rnav2(i)=rnavt2(i)
c
c  Set flag indicating that nav data is ready
c
         navl=.true.
c
c  If pointing record is not in range keep on looping
c
 60   if(.not.(rhpt1.lt.rhb .and. rhb.le.rhpt2)) goto 80
c
c  If pointing record is in range store data in final variables
c
         rhp1=rhpt1
         rhp2=rhpt2
         do 70 i=1,9
            rpt1(i)=rptt1(i)
 70         rpt2(i)=rptt2(i)
c
c  Set flag indicating that pointing data is ready
c
         ptl=.true.
c
c  Read next record if either nav or pnt records are out of range
c
 80   if(.not.navl .or. .not.ptl) goto 20
c
c  Interpolate nav and pnt data to the requested time
c  the routine INTRP is currently located in the module phase2.f
c
      do 90 i=1,56
 90      call intrp(rhn1,rhn2,rhb,rnav1(i),rnav2(i),rnav(i))
      do 100 i=1,9
 100     call intrp(rhp1,rhp2,rhb,rpt1(i),rpt2(i),rpt(i))
c
c  The final nav and pnt data is in the variables RNAV and RPT
c
c  Fill sblk - use xcru to calculate cruise relevent quantities
c
      call xcru(rnav,rpt)
c
c  Reset data ready flags
c
      navl=.false.
      ptl=.false.
      return
      end
c
c========================================================
c     END of EXECUTION
c========================================================
c
      subroutine xcru(rnav,rpt)
c
c  Compute trajectory parameters and fill the common block
c  SBLK and elements in ROTMAX and ANSWER
c
c  MODIFIED from xspn(rnav,rpt,rtb,iplnt)   in phase2.f file
c  Revised at 02/16/94  by Adam Szabo
c
c  Dimension input variables:
c
      dimension rnav(56), rpt(9)
c
c  RNAV: Navigation data interpolated beween SEDR points above.
c   1-6:   Cartesian state of S/C in Earth centered ECL50 (km, km/s)
c   7-12:  Cartesian state of S/C in Sun centered ECL50 (km, km/s)
c   13-18: Cartesian state of S/C in Jupiter centered ECL50 only
c          for times covered by Launch and Cruise data, for extended
c          Cruise it contains the S/C state in Uranus centered ECL50.
c          The two modes can be distinguised by looking at elements
c          45-56, which are filled with 0.0 for Launch and Cruise modes.
c   19-24: Cartesian state of S/C in Saturn centered ECL50 for
c          Launch and Cruise times; Neptune centered ECL50 S/C state
c          for Extended Cruise mode (km, km/s)
c   25-30: Cartesian state of Earth in Sun centered ECL50 (km, km/s)
c   31-36: Cartesian state of Jupiter in Sun centered ECL50 (km, km/s)
c   37-42: Cartesian state of Saturn in Sun centered ECL50 (km, km/s)
c   43:    Sun - S/C Range (km)
c   44:    Earth-Sun-S/C angle (degrees)
c   45-50: Cartesian state of Uranus in Sun centered ECL50 (km, km/s)
c          Only available for Extended Cruise mode. In other modes
c          it is filled with 0's.
c   51-56: Cartesian state of Neptune in Sun centered ECL50 (km, km/s)
c          Only available for Extended Cruise mode. In other modes 
c          it is filled with 0's.
c
c  RPT: Pointing data interpolated between SEDR points above.
c   1-3:  Cartesian unit vector of the S/C X-axis in ECL50
c   4-6:  Cartesian unit vector of the S/C Y-axis in ECL50
c   7-9:  Cartesian unit vector of the S/C Z-axis in ECL50
c   Taken as a 3x3 matrix, it is also the rotation matrix from
c   S/C to ECL50 coordinates.
c
c  Dimension internal variables:
c  S/C position and velocity in Sun centered ECL50
c
      dimension srecl(3),svecl(3)
c
c  Position and velocity vectors of planets in ECL50, sun centered.
c  in the order of Jupter, Saturn, Uranus and Neptune
c
      real*4 erecl(3),evecl(3)
      real*4 jrecl(3),jvecl(3)
      real*4 sarecl(3),savecl(3)
      real*4 urecl(3),uvecl(3)
      real*4 nrecl(3),nvecl(3)
c
c  Nominal solar wind velocity in RTN, and in aberated S/C coordinates
c
      real*4 swnom(3), swsc(3)
c
c  Component of nominal Solar Wind flow into each cup
c
      real*4 vcomp(4)
c
c  S/C velocity in S/C coordinates and its reverse
c
      real*4 svsc(3),svm(3)
c
c  S/C and Earth positions in heliographic coordinates
c
      real*4 srhe(3), erhe(3)
c
c  Unit vector from S/C to Sun in ECL50 and S/C coordinates
c
      real*4 sunec(3), sunsc(3)
c
c  Angles between the cup normals and Sun direction
c
      real*4 sang(4)
c
c  S/C coordinates of the cup normals; null vector; cup responses
c
      real*4 cupn(3,4),v0(3),resp(4)
c
c  Rotation matrix from ECL50 to heliographic; and from there to RTN
c
      real*4 teclhe(9), thertn(9)
c
c  Rotation matrix from ECL50 to RTN
c
      real*4 tecrtn(9)
c
c  Rotation matrix from S/C to ECL50; from RTN to S/C
c
      real*4 tscecl(9), trtnsc(9)  
c
c  Generic X and Z axis, temporary vector
c
      real*4 xaxis(3), zaxis(3), temp(3)
c
c  Heliospheric nose direction in ECL50 and RTN coordinates
c
      real*4 vnose(3), vnosrtn(3)
c
c  Interstellar coordinate X and Z axis in RTN and S/C coordinates
c
      real*4 xnosrtn(3), znosrtn(3), xnossc(3), znossc(3)
c
c  Trajectory and coordinate transformation information common block
c
      common/spec/sblk(60)
c
c  The following elements of this common block are filled:
c   SBLK(1): Distance from the Sun to S/C (AU)
c        2 : Earth-Sun-S/C angle (degrees)
c        3 : Inertial heliographic latitude of S/C (degrees)
c        4 : Inertial heliographic longitude of S/C (degrees)
c        5 : Inertial heliographic range of S/C (AU) same as SBLK(1)
c        6 : Inertial heliographic latitude of Earth/IMP 8 (deg)
c        7 : Inertial heliographic longitude of Earth/IMP 8 (deg)
c        8 : Inertial heliographic range of Earth/IMP 8 (AU)
c     9-12 : Normalized response of all cups to nominal solar wind
c            (SBLK(17)) with no modulator or suppressor voltage.
c    13-16 : Component of the nominal solar wind velocity (SBLK(17))
c            along each cup normal (km/s)
c       17 : Nominal solar wind speed (if not available from mission
c            tape it is set to 400 km/s away from the Sun)
c       19 : Distance of S/C from solar rotation axis (AU)
c       20 : Distance of Earth from solar rotation axis (AU)
c    21-23 : Unit vector from S/C to Sun, S/C coordinates
c       24 : Distance of S/C from solar equatorial plane (AU)
c       25 : Distance of Earth from solar equatorial plane (AU)
c    26-29 : Angle between each cup normal and the sunward direction (deg)
c    40-42 : S/C velocity in Sun centered S/C coordinates (km/s)
c    44-52 : Rotation matrix from S/C to RTN coordinates
c
c  Rotation matrix common block
c
      common/rotmax/tcyssc(9),tscrtn(9),icoord,trot(9),voff(3),
     &  tnossc(9),teclsc(9)
c
c  The following elements of this common block are filled:
c   TCYSSC(1-9) : Rotation matrix from RTN to S/C coordinates
c   TSCRTN(1-9) : Rotation matrix from S/C to RTN coordinates
c   TNOSSC(1-9) : Rotation matrix from Interstelar to S/C coordinates
c   TECLSC(1-9) : Rotation matrix from ECL50 to S/C coordinates
c
c  Put in ANSWER common block to get nominal 
c  solar wind velocity from mission tape input
c
      COMMON/ANSWER/ANS(150),CURRNT(512)
c
c  The following elements of this common block are filled:
c   ANS(5-7) : Magnetic field vector from mission tape in
c              S/C coordinates (gamma)
c         8  : Magnitude of above magnetic field vector (gamma)
c
c  Plasma data input type common block
c  It is made available for redundant type checking for mission tapes
c
      integer*4 ntype
      character*4 jtype
      common/tptype/jtype,ntype
c
c  Initialize null vector
c
      data v0/0.,0.,0./
c
c  The cup normal vectors in S/C coordinates pointing away from cups
c   CUPN(1-3,1)=CUPN(1-3): A cup normal
c   CUPN(1-3,2)=CUPN(4-6): B cup normal
c   CUPN(1-3,3)=CUPN(7-9): C cup normal
c   CUPN(1-3,4)=CUPN(10-12): D cup normal
c
      data cupn/-.2962, -.1710, -.9397, .2962, -.1710, -.9396,
     *          0.0, .3420, -.9397, .7310, .6816, -.0349/
c
c  Definition of 1 AU in km - see comment below
c
      data AU/1.495979E8/
c
c  Conversion factor to go from degrees to radians
c
      data rd/.0174533/
c
c  Longitude of ascending node and tilt angle of sun 
c  in ECL50 (See Voyager Memo #33 for more info)
c
      DATA PHIS/75.07/, TS/7.25/
c
c  Nominal solar wind velocity vector in RTN to be used if not
c  available from mission tape (km/s)
c
      data swnom/400.,0.,0./
c
c  Set ECL50 direction of the heliosphere nose for heliosheath
c  coordinate system. This data is derived from direction
c  of the interstellar wind as reported by R. Lallement, J. L.
c  Bertaux, E. Chassefiere, and B. R. Sandel in "Lyman-Alpha
c  observations from Voyager (1-18 AU)", in Physics of the
c  Outer Heliosphere, edited by S. Grzedzielski and D.E. Page,
c  Cospar Colloquia Series, Pergamon Press, 1990. It is also
c  described in Voyager Memo #187.
c  The nose ECL50 longitude: 252 deg., latitude: 7 deg.
c    X = cos(lat) * cos(long)
c    Y = cos(lat) * sin(long)
c    Z = sin(lat)
c
      data vnose/-0.307,-0.944,0.122/
c
c  Unit axis vectors to help transformation matrix calculations.
c
      data xaxis/1.,0.,0./
      data zaxis/0.,0.,1./
c
c===================================================================
c BEGIN EXECUTION
c===================================================================
c
c
c  If Mission Tape, use true solar wind velocity if not fill value
c
      if ((ntype .eq. 4) .and. (ans(99)+ans(100)+ans(101) .lt. 1.e6))
     &   call veccpy(swnom,ans(99))
c
c  Calculate various planetary positions and velocities.
c  They are all given in units of km and km/s.
c  Loop through the three components.
c
      do 1 i=1,3
c
c  S/C Cartesion position and velocity in ECL50 (Sun centered).
c  [Position from SEDR navigation block words 13-15;
c  velocity from SEDR navigation block words 16-18].
c
         srecl(i)=rnav(i+6)
         svecl(i)=rnav(i+9)
c
c  Earth Cartesian position and velocity in ECL50 (Sun centered).
c  [Position from SEDR navigation block words 31-33;
c  velocity from SEDR navigation block words 34-36].
c
         erecl(i)=rnav(i+24)
         evecl(i)=rnav(i+27)
c
c  Jupiter Cartesian position and velocity in ECL50 (Sun centered).
c  [Position from SEDR navigation block words 37-39, launch or cruise
c  or from SEDR navigation block words 106-108, extended cruise;
c  velocity from SEDR navigation block words 40-42, launch or cruise
c  or from SEDR navigation block words 109-111, extended cruise].
c
         jrecl(i)=rnav(i+30)
         jvecl(i)=rnav(i+33)
c 
c  Saturn Cartesian position and velocity in ECL50 (Sun centered). 
c  [Position from SEDR navigation block words 43-45, launch or cruise
c  or from SEDR navigation block words 112-114, extended cruise; 
c  velocity from SEDR navigation block words 46-48, launch or cruise
c  or from SEDR navigation block words 115-117, extended cruise]. 
c 
         sarecl(i)=rnav(i+36) 
         savecl(i)=rnav(i+39)
c 
c  Uranus Cartesian position and velocity in ECL50 (Sun centered). 
c  [Position is fill of 0.0, launch or cruise (Not available)
c  or from SEDR navigation block words 37-39, extended cruise; 
c  velocity is fill of 0.0, launch or cruise (Not available)
c  or from SEDR navigation block words 40-42, extended cruise]. 
c 
         urecl(i)=rnav(i+44) 
         uvecl(i)=rnav(i+47) 
c  
c  Neptune Cartesian position and velocity in ECL50 (Sun centered).  
c  [Position is fill of 0.0, launch or cruise (Not available) 
c  or from SEDR navigation block words 43-45, extended cruise;  
c  velocity is fill of 0.0, launch or cruise (Not available)
c  or from SEDR navigation block words 46-48, extended cruise].  
c  
         nrecl(i)=rnav(i+50)  
         nvecl(i)=rnav(i+53)
c
c==================================================================
c  NO INFORMATION giving spacecraft state with respect to 
c  Jupiter or Saturn on these SEDRs, BUT this information
c  can be computed from the spacecraft and Jupiter and Saturn
c  states with respect to the sun.
c==================================================================
c
    1 continue
c
c  The input variable rpt is the rotation matrix from S/C to
c  ECL50 coordinates at the time specified by rtb. We will
c  save this matrix in TSCECL for later RTN calculations.
c
      do 2 i=1,9
 2       tscecl(i)=rpt(i)
c
c  Invert the above matrix to obtain rotation matrix from ECL50 to
c  S/C coordinates. Since rotation matrices are orthonormal (that is
c  the columns and rows are orthogonal and normal) their inverse is
c  just their transpose.
c
      CALL TEQATR(TECLSC,TSCECL)
c
c  Distance of spacecraft from sun in astronomical units (AU);
c  1 AU defined as 1.495979e8 km (ref. Allen, 3rd ed., p.16, 1973).
c  [From SEDR navigation block word 51].
c
      dscsn=rnav(43)/au
c
c  Angle between Earth-Sun vector and Spacecraft-Sun vector (degrees)
c  [From SEDR navigation block word 56].
c
      aesnsc=rnav(44)
c
c  Get unit vector from S/C to Sun in ECL50
c
      DO 556 I=1,3
 556     SUNEC(I)=-srecl(i)
      CALL VECNOR(SUNEC) 
c
c  Transform SUNEC to get unit vector to Sun in S/C coordinates
c
      call vtrns(SUNSC,teclsc,SUNEC,v0)
c
c  Find projections of this vector along cup normals to find angles
c  between cup normals and the sun in degrees.
c
      do 671 I=1,4
         SANG(I)=SUNSC(1)*CUPN(1,I)+SUNSC(2)*CUPN(2,I)+
     &    SUNSC(3)*CUPN(3,I)
         SANG(I)=ACOS(SANG(I))/RD
 671  continue
c
c  S/C velocity in S/C coordinates
c
      call vtrns(svsc,teclsc,svecl,v0)
c
c  Construct rotation matrix from S/C to heliographic  and
c  then to RTN coordinates.
c
c  From xspnur.f originally written for Uranus encounter
c  by Richard Selesnick (today is 7/25/89)
c
c  First step: Calculate transformation matrix from ECL50 to
c              to heliographic coordinates.
c
      P=PHIS*RD
      T=TS*RD
      TECLHE(1)=COS(P)
      TECLHE(2)=-SIN(P)*COS(T)
      TECLHE(3)=SIN(P)*SIN(T)
      TECLHE(4)=SIN(P)
      TECLHE(5)=COS(P)*COS(T)
      TECLHE(6)=-COS(P)*SIN(T)
      TECLHE(7)=0.
      TECLHE(8)=SIN(T)
      TECLHE(9)=COS(T)
c
c  End calculation of TECLHE: transformation matrix from
c  ECL50 to Heliographic coordinates (for DEFINITION of 
c  Heliographic coordinates, see Voyager Memorandum 33,
c  by A.J. Lazarus of Feb. 15, 1978).  This matrix 
c  is also known as M5 (Lepping's notation); note
c  that both of these systems are inertial coordinate 
c  systems.
c
c  Heliographic coordinates are also known as solar 
c  equatorial coordinates; the tilt angle between the 
c  solar rotation axis and the normal to the plane of the 
c  ecliptic is assumed to be 7.25 degrees, the longitude 
c  of the ascending node of the solar equatorial plane
c  with respect to the plane of the ecliptic is 
c  assumed to be 75.07 degrees.
c
c  N.B. In ECL50 coordinates 0 degrees longitude is 
c  the first point of Aries.  The figure in Voyager Memorandum
c  shows 0 degrees to be the location of the Earth on September 23,
c  90 degrees on December 21, 180 degrees on March 21, and 
c  270 degrees on June 22.  In the heliographic system 0 degrees
c  longitude (the X-axis) is defined to be at the ascending node
c  in this calculation.  This is consistent with the nominal
c  definition of heliographic longitude whose 0 is defined as
c  the solar meridian which passed through the ascending node
c  of the solar equator on the ecliptic on January 1, 1854
c  at Greenwich mean noon J.D. 239 8220.0; Carrington's
c  zero meridian passed the ascending node twelve hours
c  earlier (taken from p.307 of the Explanatory Supplement
c  to the Astonomical Ephemeris and The American Ephemeris 
c  and Nautical Almanac, 1961).  HOWEVER, note that the 
c  systems referred to here as Heliographic are NON-ROTATING.
c
c  [A conversion routine to get Carrington longitude
c  is in /plas0/d1/vgr/src/ibm/lib.jds/ydoyhm.f as
c  of 5/18/90].
c
c  SRECL is the S/C position in ECL50, sun centered.
c  SRHE  is the S/C position in heliographic coordinates,
c    sun centered.
c
      CALL VTRNS(SRHE,TECLHE,srecl,v0)
c
c  Calculate cylindrical radius of the S/C position from
c  the Sun's rotation axis.
c
      srhrho=sqrt(srhe(1)*srhe(1)+srhe(2)*srhe(2))
c
c  Find the Inertial Heliographic latitude and longitude
c  of the Spacecraft in degrees. Zero latitude is at the equator.
c
      call VXYZSP(SRHE,shlat,shlong,shrang)
      shlat=90.-shlat
c
c  ERECL is Earth's position in ECL50, sun centered.
c  ERHE is Earth's position in heliographic coordinates,
c   sun centered.
c
      CALL VTRNS(erhe,TECLHE,erecl,v0)
c
c  Calculate cylindrical radius of the Earth's position from
c  the Sun's rotation axis.
c
      erhrho=sqrt(erhe(1)*erhe(1)+erhe(2)*erhe(2))
c
c  Find the Inertial Heliographic latitude and longitude 
c  of Earth in degrees. Zero latitude at the equator.
c
      call VXYZSP(erhe,ehlat,ehlong,ehrang)
      ehlat=90.-ehlat
c
c  Calculate transformation matrix from inertial heliographic to
c  RTN coordinates based on S/C location.
c
      A=SQRT(SRHE(1)**2+SRHE(2)**2)
      B=SQRT(A**2+SRHE(3)**2)
      THERTN(1)=SRHE(1)/B
      THERTN(2)=-SRHE(2)/A
      THERTN(3)=-SRHE(1)*SRHE(3)/(A*B)
      THERTN(4)=SRHE(2)/B
      THERTN(5)=SRHE(1)/A
      THERTN(6)=-SRHE(2)*SRHE(3)/(A*B)
      THERTN(7)=SRHE(3)/B
      THERTN(8)=0.
      THERTN(9)=A/B
c
c  End calculation of THERTN: transformation matrix from
c  [inertial] Heliographic coordinates to [non-inertial]
c  RTN coordinates (for DEFINITION of RTN coordinates,
c  see Voyager Memorandum 33, A.J. Lazarus of Feb. 15, 1978).
c  This matrix is also known as M(theta,beta).
c 
c  Multiply the above two matrices together to obtain
c  transformation matrix from ECL50 to RTN
c
      CALL TEQAB(TECRTN,THERTN,TECLHE)
c
c  Again by multiplying obtain transformation matrix from
c  S/C to RTN coordinates.
c
      CALL TEQAB(TSCRTN,TECRTN,TSCECL)
c
c ----------------------------------------------------------------
c
c  Nominal cup response to solar wind velocity or to a 400 km/s 
c  wind blowing radially outward from the sun if real solar wind
c  data is not available from mission tape.
c
c  First get TRTNSC: rotation matrix from RTN to S/C coordinates
c                    by inverting the S/C to RTN matrix.
c
      call teqatr(trtnsc,tscrtn)
c
c  Second get SVM: the negative of the S/C velocity in S/C coordinates
c
      call vsum3(svm,-1.,svsc)
c
c  Third get SWSC: the nominal velocity vector in S/C coordinates
c                  correctly aberrated for S/C motion.
c
      call vtrns(swsc,trtnsc,swnom,svm)
c
c  Fourth get RESP: the response of each cup to this velocity
c                   assuming no modulator or suppressor voltage
c                   and normalized to normal incidence.
c
      call cuprsp(swsc,resp)
c
c  Fifth get VCOMP: the component of this velocity into each cup
c
      do 10 i=1,4
 10      vcomp(i)=-1.*vecscl(swsc,cupn(1,i))
c
c  Calculate transformation matrix from interstellar coordinate system
c  to S/C coordinates. The interstellar coordinate system is defined
c  as: X axis is the same as the RTN R axis, Z axis is |R x Vnose|
c  where Vnose is the direction of the heliosphere nose (opposite of
c  the interstellar flow vector), and Y completes the right handed
c  system.
c
c  First: Convert Vnose to RTN coordinates.
c
      call vtrns(vnosrtn,tecrtn,vnose,v0)
c
c  Second: Calculate X and Z axis in RTN coordinates:
c
      call veccpy(xnosrtn,xaxis)
      call veccrs(znosrtn,xaxis,vnosrtn)
      call vecnor(znosrtn)
c
c  Third: Express the Interstellar X and Z axis in S/C coordinates.
c
      call vtrns(xnossc,trtnsc,xnosrtn,v0)
      call vtrns(znossc,trtnsc,znosrtn,v0)
c
c  4th: Get rotation matrix which will take the generic XAXIS and
c  ZAXIS into the S/C coordinate vectors XNOSSC and ZNOSSC.
c
      call rotmtx(xaxis,zaxis,xnossc,znossc,tnossc,temp,angle)
c
c  If using a mission tape input, get magnetic field in payload (S/C)
c  coordinates
c
      if (ntype .eq. 4) call vtrns(ans(5),trtnsc,ans(96),v0)
c
c  Get magnitude of the magnetic field
c
      if (ntype .eq. 4) ans(8)=sqrt(ans(5)*ans(5)+
     +ans(6)*ans(6)+ans(7)*ans(7))
c
c------------------------------------------------------------
c
c  Fill SBLK common block
c
c  SBLK(1)      Distance from the Sun to S/C in AU
c
      sblk(1)=dscsn
c
c  SBLK(2)      Angle between Earth-Sun vector and 
c              Spacecraft-Sun vector (degrees)
      sblk(2)=aesnsc
c
c  SBLK(3)      Inertial Heliographic Latitude of Spacecraft (deg)
c
      sblk(3)=shlat
c
c  SBLK(4)      Inertial Heliographic Longitude of Spacecraft (deg)
c
      sblk(4)=shlong
c
c  SBLK(5)      Inertial Heliographic Range of Spacecraft (AU)
c               [Same as SBLK(1)]
c
      sblk(5)=shrang/AU
c
c  SBLK(6)      Inertial Heliographic Latitude of Earth/IMP 8 (deg)
c
      sblk(6)=ehlat 
c
c  SBLK(7)      Inertial Heliographic Longitude of Earth/IMP 8 (deg)
c
      sblk(7)=ehlong 
c
c  SBLK(8)      Inertial Heliographic Range of Earth/IMP 8(AU) 
c
      sblk(8)=ehrang/AU
c
c  SBLK(9-12)   Response of all cups to nominal solar wind velocity.
c  SBLK(13-16)  Component of this velocity into each cup. (km/s)
c
      do 200 i=1,4
         sblk(i+8)=resp(i)
 200     sblk(i+12)=vcomp(i)
c
c  SBLK(17)     Nominal solar wind speed (400 km/s if mission tape
c               does not have this info)
c
      sblk(17)=sqrt(swnom(1)*swnom(1)+
     &         swnom(2)*swnom(2)+swnom(3)*swnom(3))
c
c  SBLK(19)     Distance of S/C from solar rotation axis (AU)
c
      sblk(19)=srhrho/AU
c
c  SBLK(20)     Distance of Earth from solar rotation axis (AU)
c
      sblk(20)=erhrho/AU
c
c  SBLK(21-23)  Unit vector from S/C to sun, Spacecraft coordinates
c
      call veccpy(sblk(21),SUNSC)
c
c  SBLK(24)     Distance of S/C from solar equatorial plane (AU)
c
      sblk(24)=srhe(3)/AU
c
c  SBLK(25)     Distance of Earth from solar equatorial plane (AU) 
c
      sblk(25)=erhe(3)/AU
c
c  SBLK(26-29)  Angle between cup normal and direction to sun for 
c               each cup (degrees)
c
      do 300 i=1,4
 300     sblk(i+25)=sang(i)
c
c SBLK(40-42)  Spacecraft velocity in S/C coordinates, sun centered.
c              [This is in an inertial frame, planet centered
c              for phase2]. (km/s)
c
      call veccpy(sblk(40),svsc)
c
c  SBLK(44-52)  Transformation matrix from S/C to RTN coordinates
c               [In phase2 this is the transformation matrix from
c               spacecraft to cylindrical (non-rotating) coordinates
c               aligned with the planetary rotation axis].
c
      do 400 i=1,9
         sblk(i+43)=tscrtn(i)
c
c  Load TCYSSC in ROTMAX to be consistent in xanal.f
c  Should be loaded with TRTNSC
c
         tcyssc(i)=trtnsc(i)
  400 continue
c
      return
      end
c
c=============================================================
c   END EXECUTION
c=============================================================