← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/tsa/src/tsabnd.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/tsabnd.for
 
1
C @(#)tsabnd.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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.COPYRIGHT (c) 1992 European Southern Observatory & Copernicus Astron. Center
 
31
C.IDENT     tsabnd.for
 
32
C.AUTHOR    Alex Schwarzenberg-Czerny, Copernicus Astron. Center, Warsaw
 
33
C.KEYWORD   MIDAS, time series analysis, BAND/TIME
 
34
C.LANGUAGE  FORTRAN 77
 
35
C.PURPOSE   Evaluate suitable frequency band for time analysis
 
36
C.RETURNS   None
 
37
C.ENVIRON   TSA context
 
38
C.VERSION   0.0               June 1992
 
39
C 
 
40
C 021031        last modif
 
41
C 
 
42
C-----------------------------------------------------------------------------
 
43
C
 
44
C
 
45
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
46
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
47
C
 
48
      CHARACTER*60 INAME                !  NAME OF INPUT TABLE
 
49
      INTEGER      MAXOBS               !  MAXIMUM NUMBER OF SCANNED OBSERV.
 
50
C
 
51
      REAL*8       START,STEP           !  OUTPUT FREQUENCY GRID
 
52
      INTEGER      NSTEPS
 
53
C
 
54
      INTEGER      IACTS,KUN,KNUL
 
55
      INTEGER      TID,IWORK,ITIME
 
56
      INTEGER      NCOL,ICOL,NROW,IROW,ISOR
 
57
      INTEGER      LFIELD,TTYP
 
58
C 
 
59
      INTEGER*8    PTIME,PWORK
 
60
C 
 
61
      CHARACTER*10 FORM
 
62
      CHARACTER*80 TEXT
 
63
C
 
64
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
 
65
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
66
C
 
67
C
 
68
C   Get parameters
 
69
C
 
70
      CALL STSPRO ('tsabnd')
 
71
      CALL STKRDC ('IN_A',   1,1,60,IACTS,INAME ,KUN,KNUL,ISTAT)
 
72
      CALL STKRDI ('MAXOBS',   1, 1,IACTS,MAXOBS,KUN,KNUL,ISTAT)
 
73
C
 
74
C   Map input data
 
75
C
 
76
      CALL TBTOPN (INAME,F_I_MODE,TID,ISTAT)
 
77
      CALL TBIGET (TID,NCOL,NROW,ISOR,ICOL,IROW,ISTAT)
 
78
      CALL TBLSER (TID,'TIME',ITIME,ISTAT)
 
79
      IF (ITIME.LT.0) THEN
 
80
        CALL STETER(3,'Column :TIME not found')
 
81
      ENDIF
 
82
      CALL TBFGET (TID,ITIME,FORM,LFIELD,TTYP,ISTAT)
 
83
      CALL TBDGET (TID,ISTORE,ISTAT)
 
84
      IF (ISTORE.NE.F_TRANS) THEN
 
85
        TEXT='Input table '//INAME//' stored not transposed'
 
86
        CALL STETER(1,TEXT)
 
87
      ENDIF
 
88
      IF (TTYP.NE.D_R8_FORMAT) THEN
 
89
        CALL STETER(2,'Column :TIME must be in DOUBLE PRECISION')
 
90
      ENDIF
 
91
      CALL TBCMAP (TID,ITIME,PTIME,ISTAT)
 
92
C
 
93
C   Create work table
 
94
C
 
95
      CALL STFCRE('ZZMIDWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
 
96
     $    NROW,IWORK,ISTAT)
 
97
      CALL STFMAP(IWORK,F_X_MODE,1,NROW,IACTS,PWORK,ISTAT)
 
98
C
 
99
C   Evaluate frequency band
 
100
C
 
101
      CALL TIMBAND(MADRID(PTIME),MADRID(PWORK),NROW,MAXOBS,
 
102
     $    START,STEP,NSTEPS)
 
103
C
 
104
C   Return keyword values
 
105
C
 
106
      CALL STKWRD('STARTTSA', START, 1,1,KUN,ISTAT)
 
107
      CALL STKWRD('STEPTSA',  STEP,  1,1,KUN,ISTAT)
 
108
      CALL STKWRI('NSTEPS',NSTEPS,1,1,KUN,ISTAT)
 
109
      CALL STTPUT(
 
110
     $   'Keywords STARTTSA,STEPTSA and NSTEPS are set now.'
 
111
     $   ,ISTAT)
 
112
C
 
113
C   Wind-up
 
114
C
 
115
      CALL STSEPI
 
116
C
 
117
      END
 
118
C
 
119
C
 
120
C
 
121
      SUBROUTINE TIMBAND(TIM,DEL,NOBS,MAXOBS,START,STEP,NSTEPS)
 
122
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
123
      INTEGER      NOBS,NSTEPS,MAXOBS
 
124
      DOUBLE PRECISION       TIM(NOBS),DEL(NOBS)
 
125
      REAL*8       START,STEP
 
126
C
 
127
      INTEGER      NOBS1,IOBS,MAXSTEPS
 
128
      REAL*8       FINISH,STEPS,SIZE
 
129
      CHARACTER*80 TEXT
 
130
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
 
131
      DATA         MAXSTEPS/30000/
 
132
C
 
133
      IF (MAXOBS.EQ.0) MAXOBS=NOBS
 
134
      NOBS1=MIN(NOBS-1,MAXOBS)
 
135
      IF (NOBS1.LT.5) THEN
 
136
        CALL STETER(10,'Too few observations or MAXOBS too small')
 
137
      ENDIF
 
138
C
 
139
      DO 1 IOBS=1,NOBS1
 
140
        DEL(IOBS)=TIM(IOBS+1)-TIM(IOBS)
 
141
        IF (DEL(IOBS).LT.ZERO) THEN
 
142
          CALL STETER(14,':TIME must be sorted in ascending order')
 
143
        ENDIF
 
144
 1    CONTINUE
 
145
      CALL SORT(NOBS1,DEL)
 
146
C
 
147
      STEP=TIM(NOBS)-TIM(1)
 
148
      IF (STEP.LE.ZERO) THEN
 
149
        CALL STETER(11,'Input table has wrong :TIME numbers')
 
150
      ENDIF
 
151
      STEP=0.3/STEP
 
152
C
 
153
      SIZE=LOG(ONE*NOBS*ONE)
 
154
      IOBS=NOBS1*(ONE/(6.+0.3*SIZE)+0.05)+1
 
155
      FINISH=DEL(IOBS)
 
156
      IF (FINISH.LE.ZERO) THEN
 
157
        CALL STETER(12,
 
158
     $ 'Too finely spaced observations: bin them coarsly')
 
159
      ENDIF
 
160
      FINISH=HALF/FINISH*(DEL(NOBS1/2)/FINISH)**(0.6)
 
161
C
 
162
      STEPS=FINISH/STEP
 
163
C
 
164
      CALL STTPUT(' RESULTS OF FREQUENCY BAND EVALUATION:',ISTAT)
 
165
      WRITE(TEXT,'(2(A,1PE10.1))') 'Max. Frequency: ',FINISH,
 
166
     $   '  Resolution: ',STEP
 
167
      CALL STTPUT(TEXT,ISTAT)
 
168
      WRITE(TEXT,'(A,1PE10.1)') 'No. of points:  ',STEPS
 
169
      CALL STTPUT(TEXT,ISTAT)
 
170
      IF (STEPS.GT.MAXSTEPS) THEN
 
171
        STEPS=MAXSTEPS
 
172
        TEXT='*** DANGER *** Data span too long interval'//
 
173
     $     ' for good sampling of periodogrammes.'
 
174
        CALL STTPUT(TEXT,ISTAT)
 
175
        TEXT='Analysing data split into shorter'//
 
176
     $        'intervals and taking'//
 
177
     $       ' average of periodogrammes '//
 
178
     $       'will help by reducing resolution.'
 
179
        CALL STTPUT(TEXT,ISTAT)
 
180
      ENDIF
 
181
C
 
182
      NSTEPS=STEPS
 
183
      STEP=FINISH/STEPS
 
184
      START=ZERO
 
185
      END
 
186
C
 
187
C
 
188
C
 
189