PROGRAM SOLARM
C FINDS COMPONENTS OF SOLAR MOTION USING INPUTTED COORDINATES
C AND VELOCITIES OF A GROUP OF STARS
C ENTER DATA IN UNITS OF RA, DEC, & RADIAL VELOCITY
C RA AND DEC IN DECIMAL UNITS OF DEGREES
C TAKES INPUT FILE LABELLED solarm.dat
REAL RA,DEC,VR
REAL A1,A2,A3,B1,B2,B3,C1,C2,C3,D1,D2,D3,E1,E2,E3,Q
REAL DET,DETX,DETY,DETZ,X,Y,Z,DETDX,DETDY,DETDZ,DX,DY,DZ
REAL DETEX,DETEY,DETEZ,EX,EY,EZ,F1,F2,F3
DATA A1,A2,A3,B1,B2,B3 / 0.0,0.0,0.0,0.0,0.0,0.0 /
DATA C1,C2,C3,D1,D2,D3 / 0.0,0.0,0.0,0.0,0.0,0.0 /
DATA E1,E2,E3,I / 0.0,0.0,0.0,0 /
DIMENSION A(100),D(100),V(100),P(100)
OPEN (UNIT=7, FILE='solarm.dat', STATUS='OLD')
10 READ (7,*,END=12) RA, DEC, VR
I = I + 1
A(I) = RA
D(I) = DEC
V(I) = VR
A1 = A1 + COSD(A(I))*COSD(A(I))*COSD(D(I))*COSD(D(I))
A2 = A2 + COSD(A(I))*SIND(A(I))*COSD(D(I))*COSD(D(I))
A3 = A3 + COSD(A(I))*COSD(D(I))*SIND(D(I))
B1 = A2
B2 = B2 + SIND(A(I))*SIND(A(I))*COSD(D(I))*COSD(D(I))
B3 = B3 + SIND(A(I))*COSD(D(I))*SIND(D(I))
C1 = A3
C2 = B3
C3 = C3 + SIND(D(I))*SIND(D(I))
D1 = D1 - V(I)*COSD(A(I))*COSD(D(I))
D2 = D2 - V(I)*SIND(A(I))*COSD(D(I))
D3 = D3 - V(I)*SIND(D(I))
GOTO 10
12 DET = A1*(B2*C3 - B3*C2) - A2*(B1*C3 - B3*C1) + A3*(B1*C2 - B2*C1)
DETX = D1*(B2*C3 - B3*C2) - A2*(D2*C3 - B3*D3) + A3*(D2*C2 - B2*D3)
DETY = A1*(D2*C3 - B3*D3) - D1*(B1*C3 - B3*C1) + A3*(B1*D3 - D2*C1)
DETZ = A1*(B2*D3 - D2*C2) - A2*(B1*D3 - D2*C1) + D1*(B1*C2 - B2*C1)
X = DETX/DET
Y = DETY/DET
Z = DETZ/DET
DO 13 J = 1,I
Q = Z*SIND(D(J))
P(J) = 0.0 - X*COSD(A(J))*COSD(D(J)) - Y*SIND(A(J))*COSD(D(J)) - Q
E1 = E1 - ABS(V(J) - P(J))*COSD(A(J))*COSD(D(J))
E2 = E2 - ABS(V(J) - P(J))*SIND(A(J))*COSD(D(J))
E3 = E3 - ABS(V(J) - P(J))*SIND(D(J))
13 CONTINUE
F1 = E1/REAL(I - 1)
F2 = E2/REAL(I - 1)
F3 = E3/REAL(I - 1)
14 DETDX = E1*(B2*C3 - B3*C2) - A2*(E2*C3 - B3*E3) + A3*(E2*C2 - B2*E3)
DETDY = A1*(E2*C3 - B3*E3) - E1*(B1*C3 - B3*C1) + A3*(B1*E3 - E2*C1)
DETDZ = A1*(B2*E3 - E2*C2) - A2*(B1*E3 - E2*C1) + E1*(B1*C2 - B2*C1)
DETEX = F1*(B2*C3 - B3*C2) - A2*(F2*C3 - B3*F3) + A3*(F2*C2 - B2*F3)
DETEY = A1*(F2*C3 - B3*F3) - F2*(B1*C3 - B3*C1) + A3*(B1*F3 - F2*C1)
DETEZ = A1*(B2*F3 - F2*C2) - A2*(B1*F3 - F2*C1) + F1*(B1*C2 - B2*C1)
DX = DETDX/DET
DY = DETDY/DET
DZ = DETDZ/DET
EX = DETEX/DET
EY = DETEY/DET
EZ = DETEZ/DET
WRITE (6,15) X, DX, EX
15 FORMAT (1X,'X = ',F13.5,' DX = +-',F13.5,' SE = +-',F13.5)
WRITE (6,16) Y, DY, EY
16 FORMAT (1X,'Y = ',F13.5,' DY = +-',F13.5,' SE = +-',F13.5)
WRITE (6,17) Z, DZ, EZ
17 FORMAT (1X,'Z = ',F13.5,' DZ = +-',F13.5,' SE = +-',F13.5)
STOP
END