program llspa
c time in seconds, file b1d format
C THIS program FITS Y=A+BX+CSIN W(X-PZ) TO DATA
C X IS THE INDEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
C Y IS THE DEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
C N IS THE NUMBER OF POINTS: INPUT DATA
C W IS THE "frequency" OF THE SINE FUNCTION:INPUT DATA
C R IS A FOUR LINES ARRAY WHERE:R(1)=A,R(2)=B,R(3)=C,R(4)=PZ:OUTPUT DATA
C SIG2 IS THE VARIANCE OF THE FIT:SIG2=SUM((Y OBS-Y CALC)**2)*P/(N-4):
C OUTPUT DATA
C CM IS THE 4X4 COVARIANCE MATRIX:OUTPUT DATA
C THE METHOD OF FINDING A,B,C,FI IS THROUGH "CSIN W(X-PZ)=DSINWX+ECOSWX"
C USES SUBROUTINE SYMIN
C************************************************************************
C Adapted from S. O. Kepler, this version written 17 MAY 1984
C modified 4 Feb 97 to reminded to add t1
C************************************************************************
implicit double precision(a-h,o-z)
common /x/ X(250000)
common /y/ Y(250000)
dimension R(4),A(4,4),B(4),CM(4,4),sigma(4)
DIMENSION P(20),Q(20),RR(20)
CHARACTER*20 infile
pi=dacos(-1.0d0)
pi2=2.0d0*pi
321 write(*,322)
322 format(1x,'input file is (N<250,000): ', $)
read(*,323) infile
323 format(a20)
open(1,file=infile,status='old')
c *********period is given here (sec)***************
write(*,324)
324 format(' input period = ', $)
read(*,*) period
c **************************************************
C... input file
c read in data - time (x) and intensity (y)
19 i=1
read(1,*)ndata
read(1,*)t1,y(1)
x(1)=0.0d0
do 21 nd=2,ndata
i=i+1
read(1,*)tt,y(i)
x(i)=(tt-t1)
21 continue
write(*,27)ndata
27 format(' ndata= ',i5)
28 continue
CLOSE(1)
n=ndata
31 format(' finished readind data, npoints= ',i6)
C
c************************************************************
W=pi2/period
DO 1 I=1,4
DO 1 J=1,4
A(I,J)=0.0d0
B(I)=0.0d0
R(I)=0.0d0
CM(I,J)=0.0d0
1 CONTINUE
DO 2 I=1,N
XI=X(I)
YI=Y(I)
PI=1.0d0
WX=W*XI
WX=DMOD(WX,pi2)
CI=DCOS(WX)
SI=DSIN(WX)
A(1,1)=A(1,1)+PI
A(1,2)=A(1,2)+PI*XI
A(1,3)=A(1,3)+PI*SI
A(1,4)=A(1,4)+PI*CI
A(2,2)=A(2,2)+PI*XI**2
A(2,3)=A(2,3)+PI*XI*SI
A(2,4)=A(2,4)+PI*XI*CI
A(3,3)=A(3,3)+PI*SI**2
A(3,4)=A(3,4)+PI*SI*CI
A(4,4)=A(4,4)+PI*CI**2
B(1)=B(1)+PI*YI
B(2)=B(2)+PI*YI*XI
B(3)=B(3)+PI*YI*SI
B(4)=B(4)+PI*YI*CI
2 CONTINUE
DO 3 I=2,4
DO 3 J=1,I
A(I,J)=A(J,I)
3 CONTINUE
NN=4
C CALL SYMIN(A,NN)
C SUBROUTINE SYMIN (A,NN)
C
C THIS SUBROUTINE INVERTS A SYMMETRIC MATRIX.
C ***** INPUT DATA *****
C A = A SYMMETRIC MATRIX
C NN = THE RANK OF THE MATRIX. N MUST BE LESS THAN 20.
C ***** OUTPUT DATA *****
C A = THE INVERSE OF THE INPUT MATRIX
C
C DIMENSION A(NN,NN)
C DIMENSION P(20),Q(20),R(20)
1300 FORMAT(13H0SYMIN FAILED)
ZERO = 0.0d0
ONE = 1.0d0
DO 1000 M=1,NN
1000 RR(M) = ONE
DO 200 M=1,NN
BIG = ZERO
DO 3000 L=1,NN
AB = DABS(A(L,L))
IF(AB-BIG)3000,3000,4000
4000 IF(RR(L))1400,3000,1400
1400 BIG = AB
K = L
3000 CONTINUE
IF(BIG)6000,5000,6000
5000 WRITE (*,1300)
GO TO 101
C RETURN
6000 RR(K) = ZERO
Q(K) = ONE/A(K,K)
P(K) = ONE
A(K,K) = ZERO
KM1 = K-1
IF(KM1.EQ.0) GO TO 1600
DO 7000 L=1,KM1
P(L) = A(L,K)
IF(RR(L))9000,8000,9000
8000 Q(L) = A(L,K)*Q(K)
GO TO 7000
9000 Q(L) = -A(L,K)*Q(K)
7000 A(L,K) = ZERO
1600 CONTINUE
KP1 = K+1
IF(KP1.GT.NN) GO TO 1700
DO 1500 L=KP1,NN
IF(RR(L))1200,1100,1200
1200 P(L) = A(K,L)
GO TO 10000
1100 P(L) = -A(K,L)
10000 Q(L) = (-A(K,L))*Q(K)
1500 A(K,L) = ZERO
1700 CONTINUE
DO 200 L=1,NN
DO 200 K=L,NN
200 A(L,K) = A(L,K) + P(L)*Q(K)
M = NN+1
L = NN
DO 207 K=2,NN
M = M-1
L = L-1
DO 207 J=1,L
207 A(M,J) = A(J,M)
C RETURN
C NOW THE A MATRIX IS THE INVERSE OF THE INITIAL ONE
C FINDING R(I),WHERE R(1)=A,R(2)=B,R(3)=D,R(4)=E
DO 4 I=1,4
DO 4 J=1,4
R(I)=R(I)+A(I,J)*B(J)
4 CONTINUE
C FINDING SIG2
RSIG=0
DO 5 I=1,N
WX=W*X(I)
WX=DMOD(WX,pi2)
SI=DSIN(WX)
CI=DCOS(WX)
PI=1.0d0
RSIG=RSIG+((Y(I)-R(1)-R(2)*X(I)-R(3)*SI-R(4)*CI)**2)*PI
5 CONTINUE
SIG2=RSIG/(dble(N)-4.0d0)
C FINDIND THE CORRELATION MATRIX CM
DO 6 I=1,4
DO 6 J=1,4
CM(I,J)=SIG2*A(I,J)
6 CONTINUE
C C SIN W(X-PZ)=C(SIN(WX)COS(WPZ)-COS(WX)SIN(WPZ))=D SIN(WX)+E COS(WX)
C D=C COS(WPZ) E=-CSIN(WPZ)
C DEFINING R(3)=C,R(4)=PZ
PZ=DATAN2(-1.0d0*R(4),R(3))/W
C=DSQRT(R(4)**2+R(3)**2)
c332=(r(3)**2*cm(3,3)+2*(r(3)*r(4)*cm(3,4))+
&r(4)**2*cm(4,4))/(r(3)**2+r(4)**2)
67 c33=dsqrt(c332)
cm(1,3)=c33*cm(1,3)/dsqrt(cm(3,3))
cm(2,3)=c33*cm(2,3)/dsqrt(cm(3,3))
c442=(cm(4,4)*r(3)**2-2*(r(3)*r(4)*cm(3,4))+
&cm(3,3)*r(4)**2)/((r(3)**2+r(4)**2)**2*w**2)
c44=dsqrt(c442)
cm(1,4)=c44*cm(1,4)/dsqrt(CM(4,4))
cm(2,4)=c44*cm(2,4)/dsqrt(CM(4,4))
cm(3,4)=c33*c44
cm(3,3)=c332
cm(4,4)=c442
do 3910 i=2,4
do 3910 j=1,i
cm(i,j)=cm(j,i)
3910 continue
r(3)=c
r(4)=pz
do 3915 i=1,4
sigma(i)=dsqrt(cm(i,i))
3915 continue
c NO correction of the amplitude due to tapering of data
c r(3)=r(3)/0.9d0
if(r(4).lt.0.0)r(4)=r(4)+period
c calculation of the time of maximum
pmax=r(4)+period/4.0d0
if(pmax.gt.period)pmax=pmax-period
100 format('Fitting a sine curve to the time variation'/,
&'y=a+bx+c sin(wt-pz)'/,'a = ',1pe12.5,/'b = ',1pe12.5,
&/'c = ',1pe12.5,/'pz=',0pf8.3,/'period = ',0pf10.5,
&/'sig2= ',1pe12.3)
c write header*******************************
c *******************************************
432 write(*,433)
433 format(1x,'output file name is: ', $)
read(*,323) infile
open(3,file=infile,status='unknown')
write(3,'(a20)')infile
write(3,99)N
99 format('npoints= ',i6)
c output results
write(3,100)r,period,sig2
write(3,110)(sigma(i),i=1,4),t1,pmax
110 format('siga =',1pe12.3,/'sigb =',1pe12.3,
&/'sigc =',1pe12.3,/'sigpz =',0pf9.4,/
&'time of maximum (after t1=',0pf15.4,') = ',0pf9.4)
write(3,107)((cm(i,j),j=1,4),i=1,4)
107 format(31x,'Correlation Matrix'/,
&29x,'a, b, c, time of zero'/,
&1x,4(4(0pg18.6,2x)/,1x))
101 continue
end