← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/tsa/src/tsaaov.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/tsaaov.for
 
1
C @(#)tsaaov.for        19.1 (ESO-DMD) 02/25/03 13:33:25
 
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  @(#)tsaaov.for       5.1 (ESO-IPG) 4/5/93 15:58:45
 
30
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C.COPYRIGHT (c) 1992 European Southern Observatory & Copernicus Astron. Center
 
32
C.IDENT     tsaaov.for
 
33
C.AUTHOR    Alex Schwarzenberg-Czerny, Copernicus Astron. Center, Warsaw
 
34
C.KEYWORD   MIDAS, time series, AOV/TSA
 
35
C.LANGUAGE  FORTRAN 77
 
36
C.PURPOSE   Compute analysis of variance periodogramme
 
37
C.RETURNS   None
 
38
C.ENVIRON   TIME/ANAL context
 
39
C.VERSION   0.0               June 1992
 
40
C
 
41
C 021031        last modif
 
42
C 
 
43
C-----------------------------------------------------------------------------
 
44
C
 
45
C
 
46
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
47
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
48
C
 
49
      INTEGER      NAXIS
 
50
      PARAMETER    (NAXIS=2)
 
51
      REAL*8       START(NAXIS)         !  START FREQUENCY
 
52
      REAL*8       STEP(NAXIS)          !  STEP IN FREQUENCY
 
53
      INTEGER      LENGTH(NAXIS)        !  NUMBER OF FREQUENCIES
 
54
      INTEGER      ORDER                !  NUMBER OF BINS
 
55
      INTEGER      COVER                !  NUMBER OF BIN COVERAGES
 
56
      CHARACTER*60 INAME                !  NAME OF INPUT TABLE
 
57
      CHARACTER*60 ONAME                !  NAME OF OUTPUT IMAGE
 
58
C
 
59
      INTEGER      IACTS,KUN,KNUL
 
60
      INTEGER      TID,ITIME,IVAL
 
61
      INTEGER      NCOL,ICOL,NROW,IROW,ISOR
 
62
      INTEGER      LFIELD,TTYP,VTYP
 
63
      INTEGER      IDPER
 
64
      INTEGER*8    PTIME,PVAL,PPER
 
65
C 
 
66
      CHARACTER*10 FORM
 
67
      CHARACTER*80 TEXT
 
68
C
 
69
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
 
70
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
71
C
 
72
C   Get parameters
 
73
C
 
74
      CALL STSPRO ('tsaaov')
 
75
      CALL STKRDC ('IN_A'  ,1,1,60,IACTS,INAME    ,KUN,KNUL,ISTAT)
 
76
      CALL STKRDC ('OUT_A', 1,1,60,IACTS,ONAME    ,KUN,KNUL,ISTAT)
 
77
      CALL STKRDD ('STARTTSA',
 
78
     $     1, 1,IACTS,START(1) ,KUN,KNUL,ISTAT)
 
79
      CALL STKRDD ('STEPTSA',
 
80
     $     1, 1,IACTS,STEP(1)  ,KUN,KNUL,ISTAT)
 
81
      CALL STKRDI ('NSTEPS',  1, 1,IACTS,LENGTH(1),KUN,KNUL,ISTAT)
 
82
      CALL STKRDI ('ORDERTSA',
 
83
     $     1, 1,IACTS,ORDER, KUN,KNUL,ISTAT)
 
84
      CALL STKRDI ('COVER',   1, 1,IACTS,COVER, KUN,KNUL,ISTAT)
 
85
      START(NAXIS)=  ZERO
 
86
      STEP(NAXIS)=   ONE
 
87
      LENGTH(NAXIS)= 2
 
88
C
 
89
C   Map input data
 
90
C
 
91
      CALL TBTOPN (INAME,F_I_MODE,TID,ISTAT)
 
92
      CALL TBIGET (TID,NCOL,NROW,ISOR,ICOL,IROW,ISTAT)
 
93
      CALL TBLSER (TID,'TIME',ITIME,ISTAT)
 
94
      IF (ITIME.LT.0) THEN
 
95
        CALL STETER(3,'Column :TIME not found')
 
96
      ENDIF
 
97
      CALL TBLSER (TID,'VALUE',IVAL,ISTAT)
 
98
      IF (IVAL.LT.0) THEN
 
99
        CALL STETER(3,'Column :VALUE not found')
 
100
      ENDIF
 
101
      CALL TBFGET (TID,ITIME,FORM,LFIELD,TTYP,ISTAT)
 
102
      CALL TBFGET (TID,IVAL, FORM,LFIELD,VTYP,ISTAT)
 
103
      CALL TBDGET (TID,ISTORE,ISTAT)
 
104
      IF (ISTORE.NE.F_TRANS) THEN
 
105
        TEXT='Input table '//INAME//' stored not transposed'
 
106
        CALL STETER(2,TEXT)
 
107
      ENDIF
 
108
      IF (TTYP.NE.D_R8_FORMAT.OR.VTYP.NE.D_R8_FORMAT) THEN
 
109
        CALL STETER(1,'Column(s) must be of DOUBLE PRECISION type')
 
110
      ENDIF
 
111
      CALL TBCMAP (TID,ITIME,PTIME,ISTAT)
 
112
      CALL TBCMAP (TID,IVAL, PVAL, ISTAT)
 
113
C
 
114
C   Map output data
 
115
C
 
116
      CALL STIPUT (ONAME,D_R8_FORMAT,F_IO_MODE,
 
117
     $   F_IMA_TYPE,NAXIS,LENGTH,
 
118
     $   START,STEP,ONAME,'AXIS: 1/TIME DATA: UNITLESS',
 
119
     $   PPER,IDPER,ISTAT)
 
120
C
 
121
C   Compute AOV periodogramme
 
122
C
 
123
      CALL TIMAOV (MADRID(PTIME),MADRID(PVAL),
 
124
     $    MADRID(PPER),NROW,LENGTH(1),
 
125
     $   START(1),STEP(1),ORDER,COVER)
 
126
C
 
127
C   Wind-up
 
128
C
 
129
      CALL DSCUPT(IDPER,IDPER,' ',ISTAT)
 
130
      CALL STSEPI
 
131
C
 
132
      END
 
133
C
 
134
C
 
135
C
 
136
      SUBROUTINE TIMAOV(T,X,Y,N,NFREQ,FREQIN,FSTEP,NORDER,NCOVER)
 
137
C
 
138
C     Analysis of Variance Periodogramme
 
139
C     Reference:   M.N.R.A.S. 241, 153
 
140
C     Alex Schwarzenberg-Czerny               Warsaw Observatory 1988
 
141
C
 
142
C
 
143
      INCLUDE  'MID_REL_INCL:TSA_DEF.INC'
 
144
C
 
145
      INTEGER  NORDER                   ! ORDER (NUMBER OF BINS)
 
146
      INTEGER  NCOVER                   ! NUMBER OF BIN COVERS
 
147
      INTEGER  N                        ! NUMBER OF OBSERVATIONS
 
148
      INTEGER  NFREQ                    ! NUMBER OF FREQUENCIES
 
149
      DOUBLE PRECISION   T(N),X(N)      ! TIMES AND VALUES OF OBSERVATIONS
 
150
      DOUBLE PRECISION   Y(NFREQ,2)     ! PERIODOGRAMME
 
151
      REAL*8   FREQIN,FSTEP             ! FREQUENCY START AND INCREMENT
 
152
C
 
153
      INTEGER  MAXORD
 
154
      PARAMETER(MAXORD=100)
 
155
C
 
156
      INTEGER      ICOVER,IF,I,NO,IBIN,NR
 
157
      INTEGER      NINCMPL,NINSUFF,NALL
 
158
      INTEGER      NBIN(MAXORD)
 
159
      DOUBLE PRECISION       BMEAN(MAXORD)
 
160
      DOUBLE PRECISION       T0,XM,VA,SAVE,COVER,VAR
 
161
      DOUBLE PRECISION       FREQ,XMEAN,S1,S2
 
162
      CHARACTER*80 TEXT
 
163
C
 
164
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
 
165
C
 
166
      IF (NORDER.GT.MAXORD) THEN
 
167
        WRITE(TEXT,'(A,I3)')
 
168
     1     'WARNING: Used maximum alllowed order of ',MAXORD
 
169
        CALL STTPUT(TEXT,ISTAT)
 
170
        NORDER=MAXORD
 
171
      ENDIF
 
172
C
 
173
C   Set variables (incl. mean and variance)
 
174
C
 
175
             NINCMPL=0
 
176
             NINSUFF=0
 
177
             T0=T(1)
 
178
             XM=ZERO
 
179
             VA=ZERO
 
180
             DO 10 I=1,N
 
181
               XM=XM+X(I)
 
182
 10          CONTINUE
 
183
             XM=XM/N
 
184
             DO 11 I=1,N
 
185
               SAVE=X(I)-XM
 
186
               VA=VA+SAVE*SAVE
 
187
 11          CONTINUE
 
188
             DO 12 IF=1,NFREQ
 
189
               DO 14 I=1,2
 
190
                 Y(IF,I)=ZERO
 
191
 14            CONTINUE
 
192
 12          CONTINUE
 
193
C
 
194
C   Loop over coverages and frequencies
 
195
C
 
196
             DO 1 ICOVER=1,NCOVER
 
197
               COVER=(ICOVER-1.)/NCOVER
 
198
               NO=NORDER
 
199
               DO 7 IF=1,NFREQ
 
200
                 FREQ=((IF-1)*FSTEP+FREQIN)*NO
 
201
                 DO 3 IBIN=1,NO
 
202
                   BMEAN(IBIN)=ZERO
 
203
                   NBIN(IBIN)=0
 
204
 3               CONTINUE
 
205
C
 
206
C   Assign data to bins and add
 
207
C
 
208
                 DO 2 I=1,N
 
209
                   IBIN=1+MOD(INT((T(I)-T0)*FREQ+COVER),NO)
 
210
                   BMEAN(IBIN)=X(I)+BMEAN(IBIN)
 
211
                   NBIN(IBIN)=1+NBIN(IBIN)
 
212
 2               CONTINUE
 
213
C
 
214
C   Correct for incomplete bins
 
215
C
 
216
                 NR=NO
 
217
                 NALL=N
 
218
                 XMEAN=0.
 
219
                 VAR=VA
 
220
                 DO 4 IBIN=1,NO
 
221
                   BMEAN(IBIN)=BMEAN(IBIN)-NBIN(IBIN)*XM
 
222
                   IF (NBIN(IBIN).LE.1) THEN
 
223
                     NINCMPL=NINCMPL+1
 
224
                     NR=NR-1
 
225
                     XMEAN=XMEAN-BMEAN(IBIN)
 
226
                     VAR=VAR-BMEAN(IBIN)*BMEAN(IBIN)
 
227
                     NALL=NALL-NBIN(IBIN)
 
228
                   ENDIF
 
229
 4               CONTINUE
 
230
C
 
231
C    Calculate A.O.V. statistics for a given frequency
 
232
C
 
233
                 IF (NR.GT.1) THEN
 
234
                   XMEAN=XMEAN/NALL
 
235
                   S1=ZERO
 
236
                   DO 5 IBIN=1,NO
 
237
                     IF (NBIN(IBIN).GT.1) THEN
 
238
                       SAVE=BMEAN(IBIN)/NBIN(IBIN)-XMEAN
 
239
                       S1=SAVE*SAVE*NBIN(IBIN)+S1
 
240
                     ENDIF
 
241
 5                 CONTINUE
 
242
                   IF (S1.LT.VAR) THEN
 
243
                     S2=VAR-S1
 
244
                     Y(IF,1)=Y(IF,1)+S1/(NR-1)/S2*(NALL-NR)
 
245
                   ELSE
 
246
                     NINSUFF=NINSUFF+1
 
247
                     Y(IF,1)=Y(IF,1)+ONE
 
248
                     Y(IF,2)=MAX(Y(IF,2),ONE)
 
249
                     IF (NINSUFF.LE.1) THEN
 
250
                       CALL STTPUT(
 
251
     $                 'fatal calculation/rounding error',ISTAT)
 
252
                     ENDIF
 
253
                   ENDIF
 
254
                 ELSE
 
255
                   NINSUFF=NINSUFF+1
 
256
                   Y(IF,1)=Y(IF,1)+ONE
 
257
                   Y(IF,2)=MAX(Y(IF,2),TWO)
 
258
                 ENDIF
 
259
 7             CONTINUE
 
260
 1           CONTINUE
 
261
             DO 20 IF=1,NFREQ
 
262
               Y(IF,1)=Y(IF,1)/NCOVER
 
263
 20          CONTINUE
 
264
             IF (NINCMPL.GT.0.OR.NINSUFF.GT.0) THEN
 
265
               CALL STTPUT(
 
266
     $        'Statistics of bad events (see quality row for details):',
 
267
     $         ISTAT)
 
268
               WRITE(TEXT,'(2(i8,a))')
 
269
     $           NINCMPL,' underpopulated bins'
 
270
               CALL STTPUT(TEXT,ISTAT)
 
271
               WRITE(TEXT,'(2(i8,a))')
 
272
     $           NINSUFF,' of',NFREQ*NCOVER,
 
273
     $          ' bad frequencies'
 
274
               CALL STTPUT(TEXT,ISTAT)
 
275
             ENDIF
 
276
      Y(NFREQ,2)=TWO
 
277
      END
 
278
C