← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to contrib/tsa/src/tsasin.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/tsa/src/tsasin.for
 
1
C @(#)tsasin.for        19.1 (ESO-DMD) 02/25/03 13:33:26
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C.COPYRIGHT (c) 1992 European Southern Observatory & Copernicus Astron. Center
 
30
C.IDENT     tsasin.for
 
31
C.AUTHOR    Alex Schwarzenberg-Czerny, Copernicus Astron. Center, Warsaw
 
32
C.KEYWORD   MIDAS, time series, SINEFIT/TSA
 
33
C.LANGUAGE  FORTRAN 77
 
34
C.PURPOSE   Fit sine (Fourier) series
 
35
C.VERSION   0.0               June 1992
 
36
C.RETURNS   None
 
37
C.ENVIRON   TSA context
 
38
C 
 
39
C 021205        last modif
 
40
C 
 
41
C-----------------------------------------------------------------------------
 
42
C
 
43
C
 
44
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
45
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
46
C
 
47
      INTEGER      MXORDER
 
48
      PARAMETER   (MXORDER=20)          !  DEPENDS ON LENGTH OF OUTD
 
49
      REAL*8       FREQU                !  INITIAL PERIOD VALUE
 
50
      INTEGER      NORDER               !  NUMBER OF HARMONICS
 
51
      INTEGER      ITER                 !  NUMBER OF ITERATIONS
 
52
      CHARACTER*60 INAME                !  NAME OF INPUT TABLE
 
53
      CHARACTER*60 ONAME                !  NAME OF OUTPUT TABLE
 
54
      REAL*8       COEF(MXORDER+3)      !  NORDER COEFFICIENTS
 
55
C
 
56
      INTEGER      MODE,IACTS,KUN,KNUL
 
57
      INTEGER      TID,OTID,ITIME,IVAL
 
58
      INTEGER      NCOL,ICOL,NROW,IROW,ISOR
 
59
      INTEGER      LFIELD,TTYP,VTYP
 
60
      INTEGER*8     PTIME,PVAL,POTIME,POVAL
 
61
C 
 
62
      CHARACTER*10 FORM
 
63
      CHARACTER*80 TEXT
 
64
C 
 
65
      DOUBLE PRECISION E(MXORDER+3)           ! WORK TABLES
 
66
      DOUBLE PRECISION B(MXORDER+3,MXORDER+3) ! WORK TABLES
 
67
C
 
68
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
 
69
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
70
C
 
71
      DATA MODE/1/
 
72
C
 
73
C   Get parameters
 
74
C
 
75
      CALL STSPRO ('tsasin')
 
76
      CALL STKRDC ('IN_A', 1,1,60,IACTS,INAME, KUN,KNUL,ISTAT)
 
77
      CALL STKRDC ('OUT_A',1,1,60,IACTS,ONAME, KUN,KNUL,ISTAT)
 
78
      CALL STKRDD ('STARTTSA', 1, 1,IACTS,FREQU, KUN,KNUL,ISTAT)
 
79
      CALL STKRDI ('ORDERTSA',  1, 1,IACTS,NORDER,KUN,KNUL,ISTAT)
 
80
      NORDER=2*NORDER+2
 
81
      IF (NORDER.GT.MXORDER+1) THEN
 
82
        NORDER=(MXORDER-1)/2*2+1
 
83
        WRITE(TEXT,'(A,I4)') 
 
84
     $   'Will use max. permitted ORDER of', 
 
85
     $    (MXORDER-1)/2
 
86
        CALL STTPUT(TEXT,ISTAT)
 
87
      ENDIF
 
88
      CALL STKRDI ('ITER',  1, 1, IACTS,ITER,KUN,KNUL,ISTAT)
 
89
      IF (ITER.LT.0) THEN
 
90
        CALL STKRDD ('OUTPUTD',1,NORDER-1,IACTS,COEF(1),
 
91
     $    KUN,KNUL,ISTAT)
 
92
      ENDIF
 
93
C
 
94
C   Map input/output data
 
95
C
 
96
      CALL TBTOPN (INAME,F_IO_MODE,TID,ISTAT)
 
97
      OTID=TID
 
98
      CALL TBIGET (TID,NCOL,NROW,ISOR,ICOL,IROW,ISTAT)
 
99
      CALL TBLSER (TID,'TIME',ITIME,ISTAT)
 
100
      IF (ITIME.LT.0) THEN
 
101
        CALL STETER(3,'Column :TIME not found')
 
102
      ENDIF
 
103
      CALL TBLSER (TID,'VALUE',IVAL,ISTAT)
 
104
      IF (IVAL.LT.0) THEN
 
105
        CALL STETER(4,'Column :VALUE not found')
 
106
      ENDIF
 
107
      CALL TBFGET (TID,ITIME,FORM,LFIELD,TTYP,ISTAT)
 
108
      CALL TBFGET (TID,IVAL, FORM,LFIELD,VTYP,ISTAT)
 
109
      CALL TBDGET (TID,ISTORE,ISTAT)
 
110
      IF (ISTORE.NE.F_TRANS) THEN
 
111
        TEXT='Input table '//INAME//' stored not transposed'
 
112
        CALL STETER(1,TEXT)
 
113
      ENDIF
 
114
      IF (TTYP.NE.D_R8_FORMAT.OR.VTYP.NE.D_R8_FORMAT) THEN
 
115
        CALL STETER(2,
 
116
     $    'Data column(s) must be of DOUBLE PRECISION type')
 
117
      ENDIF
 
118
      CALL TBCMAP (TID,ITIME,PTIME,ISTAT)
 
119
      POTIME=PTIME
 
120
      CALL TBCMAP (TID,IVAL, PVAL, ISTAT)
 
121
      POVAL=PVAL
 
122
C
 
123
C   Map output data
 
124
C
 
125
C M. Peron 22-02-93
 
126
C A. Schwarzenberg-Czerny 07-05-93
 
127
C
 
128
C   Perform single nonlinear fit iteration, store coeficients and residuals
 
129
C
 
130
      COEF(NORDER)=FREQU
 
131
      CALL TIMFIT(MADRID(PTIME),MADRID(PVAL),MADRID(POVAL),NROW,
 
132
     $   MODE,NORDER,ITER,MXORDER,COEF,E,B)
 
133
      FREQU=COEF(NORDER)
 
134
      CALL STKWRD ('STARTTSA',  FREQU,  1,       1,KUN,ISTAT)
 
135
      CALL STKWRD ('OUTPUTD',COEF(1),1,NORDER-1,KUN,ISTAT)
 
136
C
 
137
C   Wind-up
 
138
C
 
139
      CALL DSCUPT(OTID,OTID,' ',ISTAT)
 
140
      CALL STSEPI
 
141
C
 
142
      END
 
143
C
 
144
C
 
145
C
 
146
C
 
147
C
 
148
C
 
149
C                 Subroutine  FIT
 
150
C
 
151
C     Purpose: Given a trial period, fit a Fourier series by non-linear
 
152
C        least squares method.
 
153
C
 
154
C     Data: an ASCII file containing 3 columns (number of line(integer),
 
155
C        time, value (both real*4)) and less than maxm-2 rows.
 
156
C        A trial value of the period, an order of the series (<maxo/2-1)
 
157
C        and a name of the data file are prompted for.
 
158
C
 
159
C     Results: Correlation coefficients of the fitted parameters, their
 
160
C        values and errors are displayed on the screen after each
 
161
C        iteration. Finish iterations by typing 0 when changes in
 
162
C        the parameters and std.dev. are negligible. After the final
 
163
C        iteration values of the series are subtracted from the data and
 
164
C        output as fitsin.res file.
 
165
C
 
166
C     Method: The data are fitted by the following series:
 
167
C        f(1)+f(2)*sin(w*t)+f(3)*cos(w*t)+...+f(2*ord+1)*cos(2*ord*w*t)
 
168
C        by least squares. After the initial solution
 
169
C        the value of the frequency w is varied as well as the coefficients so
 
170
C        that due to nonlinearity of the problem several iterations are
 
171
C        required. To continue iterations type 1.
 
172
C
 
173
C     A. Schwarzenberg-Czerny                      June, 1987
 
174
C
 
175
      SUBROUTINE TIMFIT(T,X,OX,NX,MODE,NORD,NITER,MXORD,F,E,B)
 
176
C
 
177
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
178
C
 
179
 
 
180
      INTEGER NX                 !  NUMBER OF DATA
 
181
      DOUBLE PRECISION  CRACOW             !  EXTERNAL FUNCTION
 
182
      DOUBLE PRECISION  T(NX)              !  TIME VECTOR
 
183
      DOUBLE PRECISION  X(NX)              !  INPUT DATA VECTOR
 
184
      DOUBLE PRECISION  OX(NX)             !  OUTPUT DATA VECTOR
 
185
      INTEGER MODE               !  MODE OF WORK
 
186
      INTEGER NORD               !  ORDER OF THE SERIES
 
187
      INTEGER MXORD              !  TABLE SIZE-3
 
188
      INTEGER NITER              !  NUMBER OF ITERATIONS
 
189
      DOUBLE PRECISION  F(MXORD+3)         !  FIT COEFFICIENTS (INCL. FREQ)
 
190
      DOUBLE PRECISION E(MXORD+3),B(MXORD+3,MXORD+3)   ! WORK TABLES
 
191
C
 
192
      DOUBLE PRECISION  FREQ
 
193
      LOGICAL START
 
194
      CHARACTER*80 MSG
 
195
      DOUBLE PRECISION TM,XM,TMX,TMI,DFREQ,PI2,CORRMX
 
196
      DOUBLE PRECISION PHMT,DPHMT,EPMAX, FREQSTEP
 
197
      INTEGER N,NP,I,N1,J,K,ITER
 
198
C
 
199
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
 
200
C
 
201
C    get data
 
202
C
 
203
      PI2=TWO*TWO*TWO*ATAN(ONE)
 
204
      START=NITER.GT.0
 
205
      NITER=ABS(NITER)
 
206
      ITER=0
 
207
      TM=ZERO
 
208
      XM=ZERO
 
209
      TMX=T(1)
 
210
      TMI=TMX
 
211
      DO 1 I=1,NX
 
212
        TMX=MAX(TMX,T(I))
 
213
        TMI=MIN(TMI,T(I))
 
214
        TM=TM+T(I)
 
215
 1      XM=XM+X(I)
 
216
      TM=TM/NX
 
217
      XM=XM/NX
 
218
C
 
219
C    input parameters and output results of previous iteration
 
220
C
 
221
      IF (MODE.EQ.1) THEN
 
222
        IF (NORD.LT.1) CALL STETER(5,'Wrong parameter')
 
223
        F(NORD)=PI2*F(NORD)
 
224
        NP=2
 
225
        IF (START) THEN
 
226
           WRITE(MSG,'(a,1pg20.6,a,i4)')
 
227
     $         ' START: frequency=',F(NORD)/PI2,' order=',(NORD-2)/2
 
228
           CALL STTPUT(MSG,ISTAT)
 
229
           N=NORD-1
 
230
           DFREQ=ZERO
 
231
        ELSE
 
232
          N=NORD
 
233
        ENDIF
 
234
      ELSE
 
235
        CALL STETER(2,' Wrong mode')
 
236
      ENDIF
 
237
      IF (START) THEN
 
238
        DO 5 I=1,N
 
239
 5         F(I)=ZERO
 
240
      ENDIF
 
241
C
 
242
C    iterations for non-linear fit
 
243
C
 
244
 100  CONTINUE
 
245
C
 
246
C           freq - current frequency
 
247
C           dfreq - its error
 
248
C           tm - mean time of obserwation
 
249
C           phmt - phase of base sinusoid at mean time
 
250
C           dphmt - its error
 
251
C           epmax - epoch of maximum of the base sinusoid
 
252
C           The harmonics are ignored by phase calculation
 
253
C
 
254
C
 
255
C     set condition and normal equations
 
256
C
 
257
        N1=N+1
 
258
        ITER=ITER+1
 
259
        DO 60 I=1,N1
 
260
          DO 60 J=1,N1
 
261
 60         B(J,I)=ZERO
 
262
        DO 3 I=1,NX
 
263
          F(NORD+1)=T(I)-TM
 
264
          F(NORD+2)=X(I)-XM
 
265
          IF (MODE.EQ.1) THEN
 
266
            CALL FITSIN(E,N1,F,NP)
 
267
          ENDIF
 
268
C  *** prevent overwriting data before end of iterations, A.S-C, 07-05-93
 
269
          IF (ITER.GT.NITER) THEN
 
270
            OX(I)=E(N1)
 
271
          ENDIF
 
272
          DO 3 J=1,N1
 
273
            DO 3 K=1,N1
 
274
 3            B(K,J)=B(K,J)+E(J)*E(K)
 
275
C
 
276
C     solve normal equations and evaluate results
 
277
C
 
278
        IF (ITER.LE.NITER) THEN
 
279
          IF (CRACOW(B,N1,MXORD+3,NX).NE.ZERO) THEN
 
280
C
 
281
C   Update and output coefficients
 
282
C
 
283
            CORRMX=ZERO
 
284
            DO 7 J=1,N
 
285
              F(J)=F(J)+B(J,N1)
 
286
              DO 7 I=1,J-1
 
287
 7              CORRMX=MAX(CORRMX,ABS(B(I,J)/SQRT(B(J,J)*B(I,I))))
 
288
            IF (ITER.LE.1) THEN
 
289
              CALL STTPUT(
 
290
     $        '   iter.        freq.       std.dev.,     max. corr.',
 
291
     $          ISTAT)
 
292
            ELSE
 
293
              FREQSTEP=B(NORD,N1)/PI2
 
294
              IF (FREQSTEP*(TMX-TMI).GT.TWO) THEN
 
295
                CALL STTPUT(
 
296
     $            '*** Large frequency step, divergence?',ISTAT)
 
297
              ENDIF
 
298
            ENDIF
 
299
            WRITE(MSG,'(1x,i7,3(1pe15.7),'' '')')
 
300
     $         ITER,F(NORD)/PI2,SQRT(MAX(B(N1,N1),ZERO)),CORRMX
 
301
            CALL STTPUT(MSG,ISTAT)
 
302
            IF (ITER.EQ.NITER) THEN
 
303
              CALL STTPUT('   coefficients     errors',ISTAT)
 
304
              F(1)=F(1)+XM
 
305
              DO 6 J=1,N
 
306
                WRITE(MSG,'(1x,2(1pe15.7),'' '')')
 
307
     $                F(J),SQRT(MAX(B(J,J),ZERO))
 
308
 6            CALL STTPUT(MSG,ISTAT)
 
309
              F(1)=F(1)-XM
 
310
C
 
311
C   Output ephemeris for base sinusoid
 
312
C
 
313
              IF (NORD.GT.ZERO.AND.MODE.EQ.1) THEN
 
314
                FREQ=F(NORD)/PI2
 
315
                IF(.NOT.START) DFREQ=SQRT(MAX(B(N,N),ZERO))/PI2
 
316
                PHMT=ATAN2(F(3),F(2))/PI2
 
317
                DPHMT=SQRT(F(2)*F(2)*B(3,3)+F(3)*F(3)*B(2,2))/
 
318
     $             (F(2)*F(2)+F(3)*F(3))
 
319
                EPMAX=(HALF*HALF-PHMT)/FREQ+TM
 
320
                IF (PHMT.LT.ZERO) PHMT=PHMT+ONE
 
321
                MSG='   frequency        error         '//
 
322
     $            'phase at m.t.  error        '
 
323
                CALL STTPUT(MSG,ISTAT)
 
324
                WRITE(MSG,'(1x,5(1pe15.7))') FREQ,DFREQ,PHMT,DPHMT
 
325
                CALL STTPUT(MSG,ISTAT)
 
326
                CALL STTPUT(
 
327
     $          '   epoch of max.    mean time    ',ISTAT)
 
328
                WRITE(MSG,'(1x,5(1pe15.7))') EPMAX,TM
 
329
                CALL STTPUT(MSG,ISTAT)
 
330
              ENDIF
 
331
            ENDIF
 
332
C
 
333
C      next iterations with period correction
 
334
            N=NORD
 
335
            START=.FALSE.
 
336
            GOTO 100
 
337
          ELSE
 
338
            CALL STETER(4,
 
339
     $          ' Singular equations, or wrong dimensions ')
 
340
          ENDIF
 
341
        ENDIF
 
342
C        end loop 100
 
343
      IF (MODE.EQ.1) THEN
 
344
        F(NORD)=F(NORD)/PI2
 
345
      ENDIF
 
346
      END
 
347
C
 
348
C
 
349
C
 
350
      SUBROUTINE FITSIN(E,N1,P,NP)
 
351
C
 
352
C    condition equations for fit of Fourier series,
 
353
C    last term contains residuals
 
354
C
 
355
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
356
C
 
357
      INTEGER N1,NP
 
358
      DOUBLE PRECISION E(N1),P(N1-1+NP)
 
359
C
 
360
      INTEGER J,NORD,NFR,L
 
361
      REAL SJ,CJ,T,X,FR,S,C,RHS,SAVE
 
362
C
 
363
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
 
364
      NORD=(N1-2)/2
 
365
      NFR=2*(NORD+1)
 
366
      T=P(NFR+1)
 
367
      X=P(NFR+2)
 
368
      FR=P(NFR)
 
369
      S=SIN(T*FR)
 
370
      C=COS(T*FR)
 
371
      SJ=S
 
372
      CJ=C
 
373
      E(NFR)=ZERO
 
374
      E(1)=ONE
 
375
      RHS=X-P(1)
 
376
      DO 4 J=1,NORD
 
377
        L=2*J
 
378
        E(L)=SJ
 
379
        E(L+1)=CJ
 
380
        E(NFR)=E(NFR)+J*T*(P(L)*CJ-P(L+1)*SJ)
 
381
        RHS=RHS-P(L)*SJ-P(L+1)*CJ
 
382
        SAVE=CJ*C-SJ*S
 
383
        SJ=CJ*S+SJ*C
 
384
 4      CJ=SAVE
 
385
      E(N1)=RHS
 
386
      END
 
387
C
 
388
C
 
389
C
 
390
C
 
391
C
 
392
C
 
393
C
 
394
C
 
395
C