PROGRAM STATPI
C	FINDS STATISTICAL PARALLAX (UPSILON AND TAU COMPONENTS) FROM
C	SPECIFIED INPUT DATA OF RA, DEC, VRAD, VIS.MAG., AND PROPER
C	MOTION COMPONENTS IN RA AND DEC.
C	SOLAR MOTION MUST BE SPECIFIED BY AA, DD, & VSUN (RA, DEC &
C	MAGNITUDE) IN INITIALIZATION
	REAL  X,Y,Z,XX,XY,XZ,YY,YZ,ZZ,A,D,VR,V,MA,MD,L,P,Q,U,T,VSTAR
	REAL  VSUM,SIN2L,USINL,ABST,ABSVR,PIUPS,PITAU,V0,VS,MVUPS,MVTAU
	REAL  B,C,BB,CC,SUM1,SUM2,SUM3,SUM4,S1,S2,S3,S4,T1,T2,SIGUPS,SIGTAU
	REAL  SIGMUP,SIGMTA
	INTEGER  I,J
	DIMENSION  A(100),D(100),VR(100),V(100),MA(100),MD(100),VSTAR(100)
	DIMENSION  L(100),P(100),Q(100),U(100),T(100),YY(100),YZ(100),ZZ(100)
	DIMENSION  B(100),C(100),BB(100),CC(100)
	DATA  ABSV,SUM1,SUM2,SUM3,SUM4 / 0.0,0.0,0.0,0.0,0.0 /
	DATA  I,VSUM,SIN2L,USINL,ABST,ABSVR / 0,0.0,0.0,0.0,0.0,0.0 /
	DATA  AA,DD,VSUN / 265.38356622,43.11170023,20.0465491 /

	WRITE (6,10)
10	FORMAT (  ' Enter decimal RA of Solar Motion RA, xxx.xxxxx')
	READ (5,11) AA
11	FORMAT (F15.10)
	WRITE (6,12)
12	FORMAT (  ' Enter decimal DEC of Solar Motion RA, xxx.xxxxx')
	READ (5,13) DD
13	FORMAT (F15.10)
	WRITE (6,14)
14	FORMAT (  ' Enter velocity of Solar Motion, xx.xxxxx')
	READ (5,15) VSUN
15	FORMAT (F15.10)

	OPEN (UNIT=7, FILE='stat.dat', STATUS='OLD')
20	READ (7,*,END=21) X, Y, Z, XX, XY, XZ
	I = I + 1
	A(I) = X
	D(I) = Y
	VR(I) = Z
	V(I) = XX
	MA(I) = XY
	MD(I) = XZ
	YY(I) = SIND(D(I))*SIND(DD) + COSD(D(I))*COSD(DD)*COSD(A(I) - AA)
	YZ(I) = COSD(DD)*SIND(A(I) - AA)
	ZZ(I) = SIND(DD)*COSD(D(I)) - SIND(D(I))*COSD(DD)*COSD(A(I) - AA)
	L(I) = ACOS(YY(I))
	P(I) = YZ(I)/SIN(L(I))
	Q(I) = ZZ(I)/SIN(L(I))
	U(I) = MA(I)*P(I) - MD(I)*Q(I)
	T(I) = MA(I)*Q(I) + MD(I)*P(I)
	VSTAR(I) = VR(I) + VSUN*YY(I)
	ABSV = ABSV + ABS(VSTAR(I))
	VSUM = VSUM + V(I)
	SIN2L = SIN2L + SIN(L(I))*SIN(L(I))
	USINL = USINL + U(I)*SIN(L(I))
	ABST = ABST + ABS(T(I))
	ABSVR = ABSVR + ABS(VR(I) + VSUN*YY(I))
	GOTO 20
21	PIUPS = 4.74*USINL/(VSUN*SIN2L)
	PITAU = 4.74*ABST/ABSVR
	V0 = VSUM/REAL(I)
	VS = ABSV/REAL(I)
	MVUPS = V0 + 5.0*ALOG10(PIUPS) + 5.0
	MVTAU = V0 + 5.0*ALOG10(PITAU) + 5.0
	J = I
	DO 22 I = 1,J
	B(I) = (U(I)*SIN(L(I))) - (USIN2L/REAL(J))
	C(I) = (SIN(L(I))*SIN(L(I))) - (SIN2L/REAL(J))
	BB(I) = ABS(T(I)) - (ABST/REAL(J))
	CC(I) = ABS(VR(I) + VSUN*YY(I)) - (ABSVR/REAL(J))
	SUM1 = SUM1 + (B(I)*B(I))
	SUM2 = SUM2 + (C(I)*C(I))
	SUM3 = SUM3 + (BB(I)*BB(I))
	SUM4 = SUM4 + (CC(I)*CC(I))
22	CONTINUE
	S1 = (SUM1*REAL(J))/(USINL*USINL*REAL(J-1))
	S2 = (SUM2*REAL(J))/(SIN2L*SIN2L*REAL(J-1))
	S3 = (SUM3*REAL(J))/(ABST*ABST*REAL(J-1))
	S4 = (SUM4*REAL(J))/(ABSVR*ABSVR*REAL(J-1))
	T1 = SQRT(S1 + S2)
	T2 = SQRT(S3 + S4)
	SIGUPS = PIUPS*T1
	SIGTAU = PITAU*T2
	SIGMUP = 2.17147241*SIGUPS/PIUPS
	SIGMTA = 2.17147241*SIGTAU/PITAU
	WRITE (6,23) V0, VSUN, VS
23	FORMAT (1X,'Vo = ',F9.4,'    Vsun = ',F9.4,'    |Vstar| = ',F9.4)
	WRITE (6,24) PIUPS, SIGUPS
24	FORMAT (1X,'UPSILON PI = ',F9.5,' +- ',F9.5)
	WRITE (6,25) MVUPS, SIGMUP
25	FORMAT (1X,'MV(UPSILON) = ',F9.5,' +- ',F9.5)
	WRITE (6,26) PITAU, SIGTAU
26	FORMAT (1X,'TAU PI     = ',F9.5,' +- ',F9.5)
	WRITE (6,27) MVTAU, SIGMTA
27	FORMAT (1X,'MV(TAU)     = ',F9.5,' +- ',F9.5)
	STOP
	END