C HELTRJCON - ADDS HELIOGRAPHIC INERTIAL LONGITUDE AND SOLAR ECLIPTIC
C COORDINATES TO HELIOS TRAJECTORY DATA FILES IN THE COHO DATA BASE.
C ORIGINAL VERSION: 11/18/93 J. F. COOPER/HUGHES STX
C DERIVED FROM ANOTHER CODE BY R. PARTHASARATHY/HUGHES STX
C CORRECTED VERSION: 6/14/94 JFC/HUGHES STX
C PREVIOUS CALCULATION OF SOLAR ECLIPTIC COORDINATES WAS INCORRECT.
C USED HELLON INSTEAD OF HILLON FOR HIL TO SEC TRANSFORMATION.
C PROBLEM CORRECTED.
Integer YYYY
Integer*4 JULJ1
Real*8 FJUL, WHOLOT, HELLAT,HILLON, HELLON, CIRC, AU, DOYDAT,D,
* XHE,YHE,ZHE,DIFJUL,SECLAT,SECLON,SPAU,XSE,YSE,ZSE
CHARACTER*50 INFILE,OUTFIL
WRITE(6,*) ' ENTER INPUT FILE NAME '
READ(5,2) INFILE
2 FORMAT(A50)
OPEN(UNIT=10,FILE=INFILE,STATUS='OLD')
WRITE(6,*) ' ENTER OUTPUT FILE NAME'
READ(5,2) OUTFIL
OPEN(UNIT=20,FILE=OUTFIL,STATUS='NEW')
AU = 1.496D8
CIRC = 360.0D0
D = 57.29578
DO 30 N=1,10000
READ(10,10) NY,ND,NH,NCARROT,SPAU,HELLON,HELLAT,
* ESSCSANG,ICODE
10 FORMAT(1X,I2,1X,I3,1X,I2,1X,I4,1X,F6.4,1X,F7.3,1X,F6.3,1X,
* F7.3,1X,I1)
NY, ND, NH = year, day, hour
NCARROT = Carrington rotation number
SPAU = Sun-Helios distance in AU
HELLON = Heliographic longitude (degrees)
HELLAT = Heliographic latitude (degrees)
ESSCSANG = Earth-Sun-Helios angle (degrees) - not used here
ICODE = Helios data availability code - not used here
DOYDAT= ND + NH/24.0 ! Fractional day of year
YYYY= NY + 1900.0
JULJ1 = 367*YYYY - ((7*(YYYY + ((1+9)/12)))/4)
* +275*1/9 + 1 +1721014 !Julian day at noon of Jan 1 of that year
FJUL = DFLOTJ(JULJ1)-0.5D0+DOYDAT-1.0D0 !Julian at exact time.
DIFJUL = FJUL - 2415020.0
C That 241.....was Julian day of 1 Jan,1900 noon, when right
C ascension of the ascending node of solar equator was 74.367 deg
R =74.367D0+1.40D0*(DIFJUL/36525.0) !Takes care of precession of
C equinox, which causes the right ascension angle to increase.
C
C Calculate solar rotation in degrees since reference epoch.(1 Jan 1854)
WHOLOT = (360.0D0/25.38D0)*(FJUL-2398220.0D0)
RANGLE = DMOD(WHOLOT,CIRC) ! Solar rotation in degrees from 1 Jan 1854
HILLON = HELLON + RANGLE ! Heliographic inertial longitude
IF(HILLON.GT.CIRC) HILLON = HILLON - CIRC
HILLAT = HELLAT ! Heliographic inertial latitude
C
XHIL = SPAU*DCOS(HILLAT/D)*DCOS(HILLON/D)
YHIL = SPAU*DCOS(HILLAT/D)*DSIN(HILLON/D)
ZHIL = SPAU*DSIN(HILLAT/D) ! Heliographic Inertial Cartesian Coord.
C We are ready for orienting the coordinate system, so X goes to
C ascending node of solar equator, Z is parallel to solar rot-axis
C TRANSFORMATION FROM SEC TO HIL COORDINATES
C XHIL = XSE*DCOS(R/D) + YSE*DSIN(R/D)
C YHIL = -0.992*XSE*DSIN(R/D)+0.992*YSE*DCOS(R/D)+0.126*ZSE
C ZHIL = 0.126*XSE*DSIN(R/D)-0.126*YSE*DCOS(R/D)+0.992*ZSE
C
C CONVERT HIL TO SEC COORDINATES
XSE=DCOS(R/D)*XHIL -0.992*DSIN(R/D)*YHIL +0.126*DSIN(R/D)*ZHIL
YSE=DSIN(R/D)*XHIL + 0.992*DCOS(R/D)*YHIL -0.126*DCOS(R/D)*ZHIL
ZSE= 0.0 + 0.126*YHIL + 0.992*ZHIL
C See Helicoor.for code for the transformation
SECLAT = D*DASIN(ZSE/SPAU) ! Solar Ecliptic Latitude
SECLON = D*DATAN2(YSE,XSE) ! Solar Ecliptic Longitude
IF(SECLON.LT.0.0) SECLON=SECLON+360.0
C We shall also calculate Carrington longitude, but won't list
C it in the output. We need to subtract the thousands of rotations
C of the sun since way back in 1 jan 1854.
C
C Note that the output of this program is in solar ecliptic
C coordinates with respect to the Mean of Date Ecliptic.
C Thus the direction of the Equinox changes for each new time.
C
WRITE(20,20) NY,ND,NH,NCARROT,SPAU,SECLAT,SECLON,HELLAT,
* HELLON,HILLON,ESSCSANG,ICODE
20 FORMAT(1X,I2,1X,I3,1X,I2,1X,I4,1X,F6.4,1X,6F7.2,1X,I1)
30 CONTINUE
50 STOP
END