danguyenwork
8/20/2013 - 4:41 PM

## An IDL program that calculates an asteroid's orbital elements given 3 coordinates on different nights.

An IDL program that calculates an asteroid's orbital elements given 3 coordinates on different nights.

``````

;/////////////////// 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
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

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

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

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

; This FUNCTION takes an angle in decimal degrees AND normalizes

; 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

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

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
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

; 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]

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

Data[2] = t_1
Data[5] = t0
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

; 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

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\"
close,2

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

PRINT, orbitalElements

; Constants are defined
k=.01720209894d; constant

; The orbital elements are assigned more descriptive names
a = orbitalElements[0]
e = orbitalElements[1]

; 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

PRINT,'DEC',DEC
IF DEC GT 180d then (DEC=DEC-360d); converts DEC
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
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

``````