Cacher is the code snippet organizer for pro developers

;/////////////////// ORBITAL DETERMINATION / EPHEMERIS GENERATION / F AND G SERIES ///////////////////

FUNCTION DMS2DEC, array
; This FUNCTION reads in an array of three elements and converts
; it to a decimal

    ; Because only the sign of the DEGREES element of the declination
    ; shows up negative, this FUNCTION tests to see IF the elements
    ; are negative or positive

    IF array[0] GT 0 THEN BEGIN
        answer = array[0] + array[1]/60d + array[2]/3600d
    ENDIF else BEGIN
        answer = array[0] - array[1]/60d - array[2]/3600d
    ENDELSE
    RETURN, answer
END

;*********************************************************************

FUNCTION deg2dms, decimal
;This procedure takes an angle in its decimal FORm and covert it back
;to degrees, minutes, seconds

    ;The FIX FUNCTION allows the program to take the integer part of a decimal)

    ;The Hour/Degrees part of the decimal is taken
    TempHour = fix(decimal)

    ;The truncated decimal is THEN subtracted from the original deciman and
    ;THEN multiplied by 60 to get the Minutes in DECIMAL
    TempMinutes = (decimal - TempHour) * 60d

    ;Same process takes place with Seconds as with Minutes
    TempSeconds = (TempMinutes - fix(TempMinutes)) * 60d

    degrees = TempHour

    ;The decimal Minutes and Seconds are THEN truncated to obtain integral
    ;Minutes and Seconds
    minutes = abs(fix(TempMinutes))
    seconds = abs(TempSeconds)

    dms = [degrees, minutes, seconds]
    RETURN, dms


END; dec2dms

;*********************************************************************

FUNCTION RAD2DEG, angle
; This FUNCTION converts an angle from degree to radian
    angle1 = angle / !dpi * 180d
    RETURN,angle1
END;

;*********************************************************************

FUNCTION DEG2RAD, angle
; This FUNCTION takes an angle in decimal degrees AND normalizes
; it in radians

        ; The angle is converted to radians
        angle = angle/180.0d * !dpi

        ; The MOD FUNCTIONs allow the computer to take the remainder when
        ; the angle is divided by 2pi to make sure that it is not greater
        ; than 2pi
        remainder = (angle MOD (2.0d * !dpi))

        ; Tf the remainder is 0 that means the angle is a multiple of 2pi.
        ; In that case the FUNCTION will RETURN 2pi. Otherwise it will
        ; RETURN the radian value (remainder)
        IF (remainder EQ 0.0) THEN remainder = 2.0d * !dpi
        RETURN, remainder
END ; DEG2RAD

;*********************************************************************

FUNCTION NormalizeRa, RA
; This FUNCTION normalizes RA

    ; The RA is first converted to a decimal, then multiplied by 15 to
    ; be in degree and finally converted to radian
    NewRA = DMS2DEC(RA)
    NewRA = NewRA * 15d
    NewRA = DEG2RAD(NewRA)
    RETURN, NewRA
END

;*********************************************************************

FUNCTION AngleAmbiguity, Sinangle,Cosangle

    ; This FUNCTION takes the sine AND the cosine of the angle
    ; to determine what quadrant the angle lies in AND find the
    ; correct value of the angle
    angle = 0d
    angle = asin(sinangle)
    IF sinangle GT 0d AND Cosangle GT 0d THEN angle=angle ; 1st quadrant
    IF sinangle GT 0d AND Cosangle LT 0d THEN angle=-angle + !dpi ;2nd quadrant
    IF sinangle LT 0d AND Cosangle LT 0d THEN angle=-angle + !dpi ; 3rd quadrant
    IF sinangle LT 0d AND Cosangle GT 0d THEN angle=angle + 2*!dpi ; 4th quadrant
    IF sinangle EQ 0d AND Cosangle EQ -1d THEN angle= !dpi ; Special case of the angle being 180
    IF sinangle EQ -1d THEN sinout = 270d ; Special case of the angle being 270

    RETURN, angle

END

;*********************************************************************

FUNCTION GETMAGNITUDE, vector
; This FUNCTION takes in a vector AND calculate its magnitude

    ; The number of components in the vector is calculated using the size FUNCTION
    length = size(vector,/N_Elements)

    ; The square of the magnitude us initiallized to 0
    squaredMagnitude = 0d

    ; Each components in the vector will be squared AND added to the squared magnitude
    FOR count=0,length-1 DO BEGIN
        squaredMagnitude = squaredMagnitude + vector[count]^2d
    ENDFOR

    ; The real magnitude is the square root of squaredMagnitude
    magnitude = SQRT(squaredMagnitude)

    ; The FUNCTION returns the magnitude
    RETURN,magnitude
END

;*********************************************************************

FUNCTION DOTPRODUCT, vector1,vector2
; This FUNCTION takes in two vectors AND returns the dot PROduct

    ; The number of components in each vector is calculated
    length1 = size(vector1,/N_Elements)
    length2 = size(vector2,/N_Elements)
    dimension1 = size(vector1,/N_Elements)
    dimension2 = size(vector2,/N_Elements)

    ; IF the two vectors ahve a dIFferent number of components, the FUNCTION will terminate
    IF (length1 NE length2) THEN $
        PRINT, 'These vectors DO not have the same number of elements' $
    ; Otherwise, the dot PROduct will be calculated using the ## operator
    ELSE BEGIN
        dotmatrix = 0d
        DotMatrix = total(vector2 * vector1)
        RETURN,dotmatrix
    END
END


;*************************************************************************

FUNCTION CROSSPRODUCT,vector1,vector2
; This FUNCTION takes in two vectors AND determine their cross PROduct

    ; The FUNCTION would first test to see IF the two vectors have the same dimensions
    ; AND number of elements
    length1 = size(vector1,/N_Elements)
    length2 = size(vector2,/N_Elements)
    dimension1 = size(vector1,/N_Elements)
    dimension2 = size(vector2,/N_Elements)

    CrossMatrix = DBLARR(3,1)

    IF (length1 NE length2) Or (dimension1 NE dimension2) THEN BEGIN
        PRINT, 'These vectors DO not have the same number of elements'
    ENDIF ELSE BEGIN

       crossmatrix[0] = vector1[1] * vector2[2] - vector1[2] * vector2[1]
       crossmatrix[1] = -vector1[0] * vector2[2] + vector1[2] * vector2[0]
       IF length1 EQ 3 THEN BEGIN
       crossmatrix[2] = vector1[0] * vector2[1] - vector1[1] * vector2[0]
       END
       RETURN,crossmatrix
    ENDELSE
END

;*********************************************************************

PRO READData, input, RSun_1, RSun, RSun1, RSun0Dot, Data
; This FUNCTION asks the user to enter the RA, DEC, AND time of three observation

    ; All the variables are initialized to doubles to preserve precision

    t_1 = 0d & RA_10 = 0d & RA_11 = 0d & RA_12= 0d & DEC_10 = 0d & DEC_11 = 0d &
    DEC_12 = 0d & t0 = 0d & RA00 = 0d & RA01 = 0d & RA02 = 0d & DEC00 = 0d &
    DEC01 = 0d & DEC02 = 0d & t1 = 0d & RA10 = 0d & RA11 = 0d & RA12 = 0d &
    DEC10 = 0d & DEC11 = 0d & DEC12 = 0d & RSun_10 = 0d & RSun_11 = 0d & RSun_12 = 0d
    RSun00 = 0d & RSun01 = 0d & RSun02 = 0d & RSun10 = 0d & RSun11 = 0d & RSun12 = 0d
    RSun0Dot0 = 0d & RSun0Dot1 = 0d & RSun0Dot2 = 0d

    OPENR, 1, "F:\SSP\InputWithCorrection.dat\"

    ; The format of the input file is:
    ;   observation#1   RA   DEC
    ;   observation#2   RA   DEC
    ;   observation#3   RA   DEC
    ;   observation#1   Sun Vector
    ;   observation#2   Sun Vector
    ;   observation#3   Sun Vector
    ;   observation#2   Sun Velocity

    READF,1, t_1, RA_10,RA_11,RA_12, DEC_10,DEC_11,DEC_12
    READF,1,t0, RA00,RA01,RA02, DEC00,DEC01,DEC02
    READF,1, t1, RA10,RA11,RA12, DEC10,DEC11,DEC12
    READF,1, RSun_10, RSun_11, RSun_12
    READF, 1, RSun00, RSun01, RSun02
    READF,1 ,RSun10, RSun11, RSun12
    ReadF,1, RSun0Dot0, RSun0Dot1, RSun0Dot2

    ; All the RAs and DECs are assigned to appropriate arrays and normalized

    RA_1 = [Ra_10, RA_11, RA_12]
    RA0 = [Ra00, RA01, RA02]
    RA1 = [Ra10, RA11, RA12]
    DEC_1 = [DEC_10, DEC_11, DEC_12]
    DEC0 = [DEC00, DEC01, DEC02]
    DEC1 = [DEC10, DEC11, DEC12]

    radRA_1 = NormalizeRA(RA_1)
    radRA0 = normalizeRA(RA0)
    radRA1 = normalizeRA(RA1)
    radDEC_1 = deg2rad(DMS2DEC(DEC_1))
    radDEC0 = deg2rad(DMS2DEC(DEC0))
    radDEC1 = deg2rad(DMS2DEC(DEC1))

    ; The DATA array includes all the RAs, DECs, times of the three observations

    Data[0] = radRA_1
    Data[1] = radDEC_1
    Data[2] = t_1
    Data[3] = radRA0
    Data[4] = radDEC0
    Data[5] = t0
    Data[6] = radRA1
    Data[7] = radDEC1
    Data[8] = t1

    ; The three Sun Vectors and the Sun Velocity arrays are created

    RSun_1 = DBLARR(3,1)
    Rsun = DBLARR(3,1)
    RSun1 = DBLARR(3,1)
    RSun_1 = [RSun_10, RSun_11, RSun_12]
    RSun = [RSun00, RSun01, RSun02]
    RSun1 = [rSun10, RSun11, RSun12]
    RSun0Dot = [RSun0Dot0, RSun0Dot1, RSun0Dot2]

    PRINT,'Data'
    PRINT, data, FORmat = '(d)'

    Close, 1
END

;*********************************************************************

FUNCTION TimeIntervals, Data
; This PROcedure converts the time interval of observation from solar days to Gaussian days

    ; The constant of conversion is .01720209895
    k = .01720209895d

    ; The Tau array include the three observation intervals
    Tau = DBLARR(3,1)

    ; The time intervals are converted to Gaussian days

    Tau[0] = k * (data[2] - data[5]) ;Tau-1 is calculated
    Tau[1] = k * (data[8] - data[2]) ;Tau0 is calculated
    Tau[2] = k * (data[8] - data[5]) ;Tau1 is calculated

    RETURN, tau

    PRINT, 'Time'
    PRINT, tau

END

;*********************************************************************

FUNCTION Get_pHat, RA, DEC
; This FUNCTION calculates pHat given RA and DEC

    pHat = DBLARR(3,1)
    pHat[0] = cos(RA) * cos(DEC)
    pHat[1] = sin(RA) * cos(DEC)
    pHat[2] = sin(DEC)
    RETURN, pHat
END

;*********************************************************************

PRO FindpHat,pHat_1, pHat0, pHat1, data
; This PROcedure calculates the pHat of three observations

    pHat_1 = DBLARR(3,1)
    pHat0 = DBLARR(3,1)
    pHat1 = DBLARR(3,1)

    pHat_1 = Get_PHat(data[0],data[1])
    pHat0 = Get_pHat(data[3],data[4])
    pHat1 = Get_pHat(data[6],data[7])

    PRINT, 'pHat_1', pHat_1
    PRINT, 'pHat', pHat0
    PRINT, 'phat1', phat1

END

;*********************************************************************

PRO pHatDot_pHat2Dot, Tau, pHat_1, pHat0, pHat1, pHat0Dot, pHat02Dot
; This PROcedure finds pHat and pHat Double Dot

    pHat0Dot = (tau[2]^(2d) * (pHat_1 - pHat0) - tau[0]^(2D) * (pHat1 - pHat0))/(tau[0] * tau[1] * tau[2])
    pHat02Dot = -2d * (tau[2] * (pHat_1-pHat0) - tau[0] * (pHat1 - pHat0))/(tau[0] * tau[1] * tau[2])

    PRINT, 'pHat0Dot', pHat0DOt,format ='(d)'
    PRINT, 'pHat02Dot',pHat02Dot
END

;*********************************************************************

FUNCTION Find_p, pHat0, RSun, r, pHat0Dot, pHat02Dot
; This FUNCTION calculates p

    ; The magnitude of the sun vector is calculated

    ; P can be expressed in term of A AND B as A/r^3 - B,
    ; where A AND B are defined as below
    Numer = CROSSPRODUCT(pHat0,pHat0Dot)
    Numer = DOTPRODUCT(Numer,RSun)

    Deno =  CROSSPRODUCT(pHat0, pHat0Dot)
    Deno = DOTPRODUCT(Deno, pHat02Dot)

    A = numer/deno

    B = 1d

    B = (B+ 1/328900.5d) / (GETMAGNITUDE(RSun)^3) * A

    p = A/r^(3d) - B
    RETURN, p
END

;*********************************************************************

FUNCTION Find_r, p, pHat0, RSun
; This FUNCTION calculates r

    ; The magnitude of the sun vector is calculated
    Rmag = GETMAGNITUDE(RSun)

    ; r is calculated
    pVector = p * pHat0
    r = SQRT(GETMAGNITUDE(RSun)^(2d) + p^2 - (2d) * DOTPRODUCT(RSun,pVector))

    RETURN, r
END


;*********************************************************************

PRO Findrp, pHat0, RSun, pHat0dot, pHat02dot, r0, p0, p0dot
; This FUNCTION uses iteration to calculate r0 AND p0.

    ; Arbitrary values of rold and rnew are given to initiate the while loop
    rold = 3d
    rNew = 2.5d

    while abs(rold - rnew) GT .000001d DO BEGIN
        ; The value of rnew is assigned to rold FOR testing purposes
        rold = rnew
        ; The values of pnew and rnew are calculated
        pNew = find_p(pHat0, RSun,rNew, pHat0dot, pHat02dot)
        print, 'pNew', pNew
        rNew = find_r(pNew, pHat0, RSun)
        print, 'rNew', rNew
    ENDWHILE

    r0 = rNew
    p0 = pNew
    PRINT, 'r0',r0
    PRINT, 'p0',p0
END

;*********************************************************************

FUNCTION Findp0dot, RSun, pHat0, pHat02Dot, pHat0Dot, r0
; This part calculates p0Dot knowing r0 AND p0

    Rmag = GETMAGNITUDE(RSun)
    Foo1 = 1d/r0^3 - (1+1/328900.5d) / Rmag^3
    Numer = CROSSPRODUCT(pHat0,pHat02Dot)
    Numer = DOTPRODUCT(Numer,RSun)
    Deno =  CROSSPRODUCT(pHat0, pHat0Dot)
    Deno = DOTPRODUCT(Deno, pHat02Dot)
    Foo2 = numer/deno

    p0Dot = -1d/2d * Foo1 * Foo2
    PRINT, 'p0dot', p0dot
    RETURN, p0dot
END

;*********************************************************************

PRO LightCorrection, pHat, pHat0Dot
; This procedure improves the accuracy of the orbital elements by taking into account
; the light-travel time

    k = .01720209895d
    c = 299792458d / 149598000000d * 86400d
    pHat = pHat - pHat0dot * ((pHat) / c) * k
    pHat = pHat / getmagnitude(pHat)
    PRINT, 'New pHat',pHat
END

;*******************************************************

PRO Find_r_and_rdot, rVector0, rDotVector0, Rsun, RsunDot, pHat0, pHat0Dot, p0, p0dot
; This Procedure calculates r and rdot

    rVector0 = p0 * pHat0 - RSun
    rDotVector0 = p0Dot * pHat0 + p0 * pHat0Dot - RSunDot

    PRINT,'rVector0',rVector0
    PRINT, 'rDotVector0',rDotVector0
END

;*********************************************************************

FUNCTION f_series, tau, r0, rVector0, rDotVector0
; This FUNCTION calculates the f part of the f AND g series
    f = 1d - tau^2 / (2*r0^3) + tau^3 * dotProduct(rVector0,rDotVector0) / (2d * r0^5)
    RETURN, f
END

;*********************************************************************

FUNCTION g_series, tau, r0, rVector0
; This FUNCTION calculates the g part of the f AND g series
    g = tau - tau^3 / (6 * r0^3)
    RETURN, g
END

;*********************************************************************

FUNCTION CalculaterVector, tau, rVector0, rDotVector0, r0
; This FUNCTION uses the f_series AND g_series FUNCTIONs to calculate rVector at a given time
    rVector = f_series(tau,r0, rVector0,rDotVector0) * rVector0 $
           + g_series(tau, r0, rVector0) * rDotVector0
    RETURN, rVector
END

;*********************************************************************

PRO f_AND_g_series, RSun_1, RSun1, pHat0, tau,r0, rVector0, rDotVector0
; This FUNCTION calculates the rVector of the first AND third observation

    ; The rVector in the first and third observation is approximated
    rVector_1 = calculatervector(tau[0], rVector0, rDotVector0,r0)
    rVector1 = calculatervector(tau[2],rVector0,rDotVector0,r0)

    ; The rVectors are then added to the Sun Vector to approximate pVector in the first
    ; and third approximation
    pVector_1 = rVector_1 + RSun_1
    pVector1 = rVector1 + RSun1

    ; Since the pVectors are not unit vectors any more, they must be divided by their
    ; magnitude to make them unit vectors again
    CheckpHat_1 = pVector_1 / GETMAGNITUDE(pVector_1)
    PRINT, 'pHat-1',CheckpHat_1
    CheckpHat1 = pVector1 / GETMAGNITUDE(pVector1)
    PRINT, 'pHat1',CheckpHat1
END

;*********************************************************************

FUNCTION FindOrbitalElements, rVector0, rDotVector0, tau,r0,data
; This FUNCTION calculates the six orbital elements

    ; Constants are defined
    u = 1d
    k = .01720209895d
    tilt = deg2rad(23.439184d)

    ; The semi-major axis is calculated
    a = 1d/ (2d / GETMAGNITUDE(rVector0) - DOTPRODUCT(rDotVector0, rDotVector0) / u)

    ; The eccentricity is calculated
    e = SQRT(1-GETMAGNITUDE(CROSSPRODUCT(rVector0,rDotVector0))^2 / (u * a))

    ; This matrix is used to rotate the rVector0 and rDotVector0 from Equatorial to
    ; Ecliptic
    TranMatrix = DBLARR(3,3)
    TranMatrix[0,0] = 1d
    TranMatrix[1,1] = cos(tilt)
    TranMatrix[2,1] = sin(tilt)
    TranMatrix[1,2] = -sin(tilt)
    TranMatrix[2,2] = cos(tilt)
  Print, 'TranMatrix', tranmatrix

    rvector0 = transpose(rVector0 # TranMatrix)
    PRINT, 'rvector0',rvector0
    rDotVector0 = transpose(rDotVector0 # TranMatrix)
    PRINT,'rDotVector0 ',rDotVector0

    ; The cross product of rVector0 and rDotVector0 is calculated
    h = CROSSPRODUCT(rvector0,rDotVector0)
    print, 'h', h

    ; The inclination is calculated
    i = atan(SQRT(h[0]^2+h[1]^2)/h[2])

    ; The longtitude of ascending node is calculated and checked FOR nagle ambiguity
    sinBOmega = h[0] / (GETMAGNITUDE(h) * sin(i))
    cosBOmega = -h[1] / (GETMAGNITUDE(h) * sin(i))
    BOmega = AngleAmbiguity(sinBOmega,cosBOmega)

    ; The distance since the asteroid crosses the ascending node is calculated
    cosU = (rVector0[0]*cos(BOmega) + rVector0[1] * sin(BOmega))/r0
    print, 'cosU', cosU
    sinU = rVector0[2] / (r0*sin(i))
    print, 'sinU', sinU
    U = angleambiguity(sinU,cosU)
    PRINT,'U',U

    ; The true anomaly is calculated. checked FOR angle ambiguity and normalized
    numer = a*(1.0-e^2) * dotproduct(rVector0,rDotVector0)
    deno = getmagnitude(h) * r0
    sinv = numer/(deno*e)
    cosv = (a*(1.0-e^2)/r0 - 1d)/e
    v = angleambiguity(sinv,cosv)
    IF (v GT !dpi) THEN v=v-2*!dpi
    PRINT, 'v',v

    ; The argument of perihelion is calculated by subtracting v from U
    w = U-v

    ; The eccentric anomaly is calculated
    EccenAnomaly = acos(1/e*(1-r0/a))
    IF v LT 0 THEN EccenAnomaly = 2*!dpi-EccenAnomaly
    PRINT, 'E',EccenAnomaly

    ; The mean anomaly at observation is calculated
    M = EccenAnomaly - e*sin(EccenAnomaly)
	print, 'm',m
    ; The mean anomaly FOR ephemeris is calculated using the mean anomaly at observation
    n = k/(sqrt(a)^3)
    print, 'n', n
    M0 = n * (2453600.5-data[5])
    print, 'm9' ,m0
    Mephem = M + M0
    print, 'mephem', mephem
    IF Mephem GT 2*!dpi THEN MEphem = Mephem - 2*!dpi

    ; The orbital elements are assigned to an array
    OrbitalElements = DBLARR(6,1)
    OrbitalElements[0] = a
    OrbitalElements[1] = e
    OrbitalElements[2] = rad2deg(i)
    OrbitalElements[3] = rad2deg(BOmega)
    OrbitalElements[4] = rad2deg(w)
    OrbitalElements[5] = rad2deg(Mephem)

    PRINT, orbitalElements, format ='(d)';
    RETURN, orbitalElements

END

;*********************************************************************

PRO EphemerisGeneration,orbitalElements
; This FUNCTION reads in the time and sun vector from a file to generate the ephemeris FOR
; that day

    ; Variables are initialized to preserve precision
    jdays=0d &Sun1=0d& Sun2=0d&Sun3=0d

    ; The name of the input file is read in


    ; The file is open FOR reading
    openr,2, "F:\SSP\EphemGeneFile.dat\"
    readf,2,jdays,Sun1,Sun2,Sun3
    close,2

    ; The SunVector is defined
    SunVector=[[Sun1],[Sun2],[Sun3]]

    PRINT, orbitalElements

    ; Constants are defined
    k=.01720209894d; constant
    ec=deg2rad(23.441884d); axial tilt of Earth

    ; The orbital elements are assigned more descriptive names
    a = orbitalElements[0]
    e = orbitalElements[1]
    i=deg2rad(orbitalElements[2])
    bomega=deg2rad(orbitalElements[3])
    lomega=deg2rad(orbitalElements[4])

    m=deg2rad(orbitalElements[5])

    ; The mean anomaly is calculated and normalized
    m=m-(2453600.5-jdays)*(k/(a^1.5))
    m=m mod (2*!dpi)
    IF m LT 0 then (m=m+2*!dpi)

	PRINT,'m', m

    EccenAnomaly=m; defines E

    FOR count=1,100 DO BEGIN; loops E-finding portion
       EccenAnomaly=EccenAnomaly-(m-(EccenAnomaly-e*sin(EccenAnomaly)))/(e*cos(EccenAnomaly)-1); loop
    endfor


    PRINT, 'EccenAnomaly',EccenAnomaly

    ; The position vector of the asteroid is defined
    positionVector=[[a*cos(EccenAnomaly)-a*e],[a*(1-e^2d)^.5d*sin(EccenAnomaly)],[0d]]

	print, 'positionVector',positionVector

    ; The four rotation matrices are defined
    SOmegaMatrix=[[cos(lomega),-sin(lomega),0d],[sin(lomega),cos(lomega), 0d],[0d,0d,1d]]
    aye=[[1d,0d,0d],[0d,cos(i),-sin(i)],[0d,sin(i),cos(i)]]
    BOmegaMatrix=[[cos(bomega),-sin(bomega),0d],[sin(bomega),cos(bomega),0d],[0d,0d,1d]]
    EclipticMatrix=[[1d,0d,0d],[0d,cos(ec),-sin(ec)],[0d,sin(ec),cos(ec)]]



    RotatedVector1=positionVector#SOmegaMatrix; multiplies original position vector
    RotatedVector2=RotatedVector1#aye; through rotation matrices to get resultant
    RotatedVector3=RotatedVector2#BOmegaMatrix
    tiltvec=RotatedVector3#EclipticMatrix

	print, 'tiltvec',tiltvec;

    pVector=SunVector+tiltvec; finds rho

	print, 'pVecor',pVector
    pHat=pVector/getmagnitude(pVector); finds pHat

	print, 'pHat',pHat

    DEC=asin(pHat[0,2]); solves FOR DEC
    PRINT,'DEC',DEC
    sinra=pHat[0,1]/cos(DEC)
    cosra=pHat[0,0]/cos(DEC)
    print, "sinra" , sinra, "cosra", cosra
    RA=angleambiguity(sinra,cosra); solves angle ambiguity FOR RA

    DEC=rad2deg(DEC)
    PRINT,'DEC',DEC
    IF DEC GT 180d then (DEC=DEC-360d); converts DEC
    RA=rad2deg(RA)/15d
     print, 'RA',RA

    PRINT, 'RA in decimal',RA
    PRINT, 'RA in dms',deg2dms(RA)
    PRINT, 'DEC in decimal',DEC
    PRINT, 'DEC in dms',deg2dms(DEC)

END

;zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

PRO Driver
; This is the main execution block
    ; The user is asked to specify the input file

    result = file_test("F:\SSP\InputWithCorrection.dat\")
    IF result EQ 0 THEN BEGIN
       PRINT, 'file not found'
    ENDIF ELSE BEGIN
       k = .01720209895d
       Data = DBLARR(3,3)
       ; The data is read in
       READData, input, RSUn_1, RSun, RSun1, RSun0DOt, Data

       ; The time interval is calculated
       tau = TimeIntervals(data)
       print, 'tau',tau
       ; Rho Hat's at the three osbervations are calculated
       FindpHat, pHat_1, pHat0, pHat1, data

       ; p Hat Dot and p Hat Double Dot are calculated
       pHatDot_pHat2Dot, Tau, pHat_1, pHat0, pHat1, pHat0Dot, pHat02Dot

       ; r0 and rhp0 are calculated
       Findrp, pHat0, RSun, pHat0dot, pHat02dot, r0, p0, p0dot

       ; pHat is correced for light travel
       LightCorrection, phat0, pHat0Dot

       ; Other values are re-calculated using the corrected pHat
       pHatDot_pHat2Dot, Tau, pHat_1, pHat0, pHat1, pHat0Dot, pHat02Dot
       Findrp, pHat0, RSun, pHat0dot, pHat02dot, r0, p0, p0dot

       ; rhoHat Dot is calculated
       p0dot=Findp0dot(RSun, pHat0, pHat02Dot, pHat0Dot, r0)
       PRINT, 'p0dot',p0dot

       ; r and r Dot are calculated
       Find_r_and_rdot, rVector0, rDotVector0, Rsun, Rsun0Dot, pHat0, pHat0Dot, p0, p0dot

       ; The f And g series approximate pHat_1 and pHat1
       f_AND_g_series, RSun_1, RSun1, pHat0, tau,r0, rVector0, rDotVector0
       PRINT, 'pHat0',pHat0

       ; The orbital elements are calculated
       orbitalElements = FindOrbitalElements(rVector0, rDotVector0, tau,r0,data)

       ; The ephemeris generation generates the RA and DEC at a given day
       EphemerisGeneration,orbitalElements
   ENDELSE
END