From: [email-masked] (Gonzalo Hernandez)
Date: Mon, 20 Nov 89 22:15:41 PST

C  Ben:
C This program takes the ( extracted ) wind data and converts it to
C magnetic coordinates and time.  Note that data enters as fort3.
C
      PROGRAM RAYPLOT
C *** takes Ray's data in UT, converts it to MLT and puts
C *** it all together. Labels in geomagnetic latitude.
C *** fort3 Ray's data in; fort4 write results out.
C *** Statements 1 -> 5; formats 101 -> 105
      DIMENSION TIMEUT(30),UMAG(30),UN(30),VN(30)
      CHARACTER*3 ISTA
      DATA PI180 /0.017453293/
      DO 1 I = 1,25
      BETA = 11 + I
      IF (BETA .GE. 24.) BETA = BETA -24.
      TIMEUT(I) = BETA
    1 CONTINUE
      FIFTEEN = 15.0*PI180
      COS15 = COS(FIFTEEN)
      SIN15 = SIN(FIFTEEN)
    3 READ(3,104,END=5) DATE,GLAT,GLONG
  104 FORMAT(3F10.3)
      IDATE = INT(DATE)
      SLONG = GLONG
      IF(SLONG .LT. 0.0) SLONG = 360 + SLONG
      READ(3,101,END=4)(UN(I), I= 1,25 )
      READ(3,101,END=4)(VN(I), I= 1,25 )
  101 FORMAT(5E12.4)
      COLAT = ABS(-90 - GLAT)*PI180
      COSLAT = COS(COLAT)
      SINLAT = SIN(COLAT)
      PLSL  = ( 120.4 - SLONG )* PI180
      SINPLSL = SIN(PLSL)
      COSPLSL = COS(PLSL)
      GEOMLAT = COS15*COSLAT+SIN15*SINLAT*COSPLSL
      GEOMLAT = ASIN(GEOMLAT)/PI180
      BETA = COS15*COSLAT+SIN15*SINLAT*COSPLSL
      BETA = ACOS(BETA)
      GEOMGLAT = BETA/PI180
      SIN1= SIN15*SINPLSL/SIN(BETA)
      BETA = ASIN(SIN1)
      COS1 = COS(BETA)
      STEST = SLONG
      IF(SLONG .LT. 0.0) STEST = 360 + SLONG
      SIGN = 1.0
      IF((STEST .GT. 30.) .AND. (STEST .LT. 210.)) SIGN = -1.
      DO 2 I = 1,25
      CALL UTLMT(GLAT,GLONG,IDATE,IDAY,IYR,TIMEUT(I),UMAG(I))
      UNT = UN(I)*COS1 + VN(I)*SIN1
      VNT = VN(I)*COS1 - UN(I)*SIN1
      UN(I) = UNT*SIGN
      VN(I) = VNT*SIGN
    2 CONTINUE
      WRITE(ISTA,102)
  102 FORMAT(3HSPO)
      WRITE(4,103)IDATE,IDAY,IYR,ISTA,GEOMLAT,(TIMEUT(I),UMAG(I),VN(I),
     *UN(I), I = 1,25 )
  103 FORMAT(1X,3I10,A3,F10.3,/,(2F7.2,2F10.2))
      GO TO 3
    4 CONTINUE
      WRITE(*,105)
  105 FORMAT(1X, 'Abnormal termination....')
    5 CONTINUE
      END
      SUBROUTINE UTLMT(GLAT,GGLONG,IDATE,IDAY,IYR,TIMEUT,UMAG)
      GLONG = GGLONG
      IF(GLONG .LT. 0.0) GLONG = 360 + GLONG
      CALL UTDAY(IDATE,IDAY,IYR)
      IYR = IYR + 1900
      UMAG = TIMEMAG(TIMEUT,IDAY,IYR,GLAT,GLONG)
      END
      SUBROUTINE UTDAY(IDATE,IDAY,IYR)
C *** IDATE AS 770606, IDAY AND IYR ARE OUTPUT AS 157 AND 77
      DIMENSION MO(12)
      DATA MO/ 0,31,59,90,120,151,181,212,243,273,304,334/
      IYR = IDATE/10000
      IMO = IDATE/100 -IYR*100
      IDAY = IDATE -(IDATE/100)*100
      LEAP = 0
      IF(IMO .GT. 12) IMO = 12
      IF(IMO .LE. 2) GO TO 1
      ITEST = IYR - 68
      IF ((ITEST/4)*4 .EQ. ITEST) LEAP = 1
    1 IDAY = MO(IMO) + IDAY + LEAP
      RETURN
      END
      FUNCTION TIMEMAG(TIMEUT,IDAY,IYR,GLAT,GLONG)
C *** MAGNETIC TIME; BASED ON COLE, 1963 FOR ECCENTRIC DIPOLE.
C *** TIMEUT IS UNIVERSAL TIME AS XX.XX
C *** IYR IN XXXX
C *** ALL LONGITUDES ARE POSITIVE EASTWARD.
C *** USES FUNCTION SOLDEC FOR SOLAR DECLINATION VALUE.
      DATA PI /3.141593/, XL1 /-0.5146 /,XL2 /0.8543/, XM1 /-.8481/,
     *XM2 /-.4947/, XN1 /-.1262/, XN2 /-.1596/, X / -.0577/, Y
     */.0320/, Z /.0184/
      TPI = 2.0*PI
      DTOR = PI/180
      RADHR = 15.*DTOR
      COLATG = GLAT*DTOR
      GLONGS = GLONG
      IF(GLONG .LT. 0.0) GLONGS = 360 + GLONG
      COLONGG = GLONGS*DTOR
      CC = PI/2.0 - COLATG
C *** H ( BELOW ) IS THE LONGITUDE OF THE SUN
      H=(24.-TIMEUT)*RADHR
      SH=SIN(H)
      CH=COS(H)
      SDEC=SOLDEC(TIMEUT,IDAY,IYR)
      SSD=SIN(SDEC)
      CSD=COS(SDEC)
      A1=SIN(CC)*COS(COLONGG)
      A2=SIN(CC)*SIN(COLONGG)
      A3=COS(CC)
      TP1=XL2*CSD*CH+XM2*CSD*SH+XN2*SSD
      TP2=XL1*CSD*CH+XM1*CSD*SH+XN1*SSD
      PHISUN=ATAN2(TP1,TP2)
      TP1=XL2*(A1-X)+XM2*(A2-Y)+XN2*(A3-Z)
      TP2=XL1*(A1-X)+XM1*(A2-Y)+XN1*(A3-Z)
      PHI0=ATAN2(TP1,TP2)
      AC=PHI0-PHISUN
      ACH=AC/RADHR
      IF(ACH .GE. 24.)ACH=ACH-24.
      IF(ACH .LT. 0.)ACH=ACH+24.
      TIMEMAG=ACH
      RETURN
      END
      FUNCTION SOLDEC(TIMEUT,IDAY,IYR)
      DATA PIODEG / 1.745329252E-2 /
      SECS = 3600.* TIMEUT
      CALL SUN(IYR,IDAY,SECS,GST,SLONG,SRASN,SDEC,ETIME)                      
      SOLDEC = SDEC*PIODEG
      RETURN
      END
      SUBROUTINE SUN(IYR,IDAY,SECS,GST,SLONG,SRASN,SDEC,ETIME)             
      DOUBLE PRECISION DJ, FDAY                                            
      DATA  RAD / 57.29578/                                                
      COTAN(A)  = 1.0 / TAN (A)                                            
      IF (IYR .LT. 1901 .OR. IYR .GT. 2099) RETURN                         
      FDAY = SECS / 86400.                                                 
      TS = FDAY*360. + 180.                                                
      IF (TS.GE. 360.) TS = TS - 360.                                      
      DJ = 365 * (IYR - 1900) + (IYR - 1901) / 4 + IDAY + FDAY - 0.5D0     
      T = DJ / 36525.                                                      
      VL = DMOD (279.696678D0+ 0.9856473354D0* DJ, 360.D0)                 
      GST=DMOD(279.690983D0+0.9856473354D0* DJ+360.D0*FDAY+ 180.D0, 360.   
     1D0)                                                                  
      G = DMOD(358.475845D0+0.985600267D0* DJ, 360.D0) / RAD               
      SLONG = VL + (1.91946 - 0.004789 * T) * SIN (G) + 0.020094 * SIN (   
     12. * G)                                                              
      OBLIQ = (23.45229 - 0.0130125 * T) / RAD                             
      SLP = (SLONG - 0.005686) / RAD                                       
      SIND = SIN (OBLIQ) * SIN (SLP)                                       
      COSD = SQRT (1. - SIND * * 2)                                        
      SDEC = RAD * ATAN (SIND / COSD)                                      
      SRASN = 180. - RAD * ATAN2 (COTAN (OBLIQ) * SIND / COSD, - COS (SL   
     1P) / COSD)                                                           
      H = GST-SRASN                                                        
      IF ( H . LT. 0.0) H = H + 360.                                       
      ETIME = (H -TS)                                                      
      RETURN                                                               
      END