← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/tsa/src/tsaint.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/tsaint.for
 
1
C @(#)tsaint.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     tsaint.for
 
32
C.AUTHOR    Alex Schwarzenberg-Czerny, Copernicus Astron. Center, Warsaw
 
33
C.KEYWORD   MIDAS, time analysis, INTERPOL/TIME
 
34
 
 
35
C.LANGUAGE  FORTRAN 77
 
36
C.PURPOSE   Interpolate an unevensampled series using its covariance function
 
37
C           Reference:  ApJ 385, 404
 
38
C.RETURNS   None
 
39
C.ENVIRON   TSA context
 
40
C.VERSION   0.0               June 1992
 
41
C 
 
42
C 021031        last modif
 
43
C 
 
44
C-----------------------------------------------------------------------------
 
45
C
 
46
C
 
47
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
 
48
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
49
C
 
50
      CHARACTER*60     INAME1         !  NAME OF 1ST OBSERVATION TABLE
 
51
      CHARACTER*60     INAME2         !  NAME OF 2ND (OUTPUT) TABLE
 
52
      CHARACTER*3      CFUNC          !  TYPE OF ACF
 
53
      DOUBLE PRECISION PARM(12)       !  PARAMETERS OF ACF
 
54
      DOUBLE PRECISION ACFLIN,ACFPOL,ACFPOW,ACFEXP,
 
55
     $                 ACFLOG,ACFIPO
 
56
      EXTERNAL         ACFLIN,ACFPOL,ACFPOW,ACFEXP,
 
57
     $                 ACFLOG,ACFIPO
 
58
      DOUBLE PRECISION TSADELUR0,TSADELUR1,TSADELUR2,TSADELUR3,
 
59
     $                 TSADELUR4,TSADELUR5,TSADELUR6,TSADELUR7,
 
60
     $                 TSADELUR8,TSADELUR9
 
61
      EXTERNAL         TSADELUR0,TSADELUR1,TSADELUR2,TSADELUR3,
 
62
     $                 TSADELUR4,TSADELUR5,TSADELUR6,TSADELUR7,
 
63
     $                 TSADELUR8,TSADELUR9
 
64
C
 
65
      INTEGER           NOBS1, NOBS2, ISIZE, ASIZE
 
66
      INTEGER           MODE,  IACTS, KUN,   KNUL
 
67
      INTEGER           TID1,  TID2,  ITIME, IDAT, IVAR
 
68
      INTEGER           NCOL,  ICOL,  NROW,  ISOR
 
69
      INTEGER           LFIELD,TTYP,  DTYP,  VTYP
 
70
      INTEGER           ICOBS, INU,   ICT
 
71
C 
 
72
      INTEGER*8         POBS1T,POBS1D,POBS1V
 
73
      INTEGER*8         POBS2T,POBS2D,POBS2V
 
74
      INTEGER*8         PCOBS,PNU,PCT
 
75
C 
 
76
      DOUBLE PRECISION  AVER,VAR
 
77
      CHARACTER*10 FORM
 
78
      CHARACTER*80 TEXT
 
79
C
 
80
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
 
81
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
82
C
 
83
C
 
84
C   Get parameters
 
85
C
 
86
      CALL STSPRO ('tsaint')
 
87
      CALL STKRDC ('IN_A',  1,1,60,IACTS,INAME1,KUN,KNUL,ISTAT)
 
88
      CALL STKRDC ('OUT_A',  1,1,60,IACTS,INAME2,KUN,KNUL,ISTAT)
 
89
      CALL STKRDC ('CFUNC',  1,1, 3,IACTS,CFUNC, KUN,KNUL,ISTAT)
 
90
      CALL STKRDD ('INPUTD',   1,12,IACTS,PARM,  KUN,KNUL,ISTAT)
 
91
C
 
92
C   Map input data
 
93
C
 
94
      CALL TBTOPN (INAME1,F_I_MODE,TID1,ISTAT)
 
95
      CALL TBIGET (TID1,NCOL,NOBS1,ISOR,ICOL,NROW,ISTAT)
 
96
      CALL TBLSER (TID1,'TIME' ,ITIME,ISTAT)
 
97
      IF (ITIME.LT.0) THEN
 
98
        CALL STETER(5,'Column :TIME not found in 1st table')
 
99
      ENDIF
 
100
      CALL TBLSER (TID1,'VALUE',IDAT ,ISTAT)
 
101
      IF (IDAT.LT.0) THEN
 
102
        CALL STETER(6,'Column :VALUE not found in 1st table')
 
103
      ENDIF
 
104
      CALL TBLSER (TID1,'VAR'  ,IVAR, ISTAT)
 
105
      IF (IVAR.LT.0) THEN
 
106
        CALL STETER(7,'Column :VAR not found in 1st table')
 
107
      ENDIF
 
108
      CALL TBFGET (TID1,ITIME,FORM,LFIELD,TTYP,ISTAT)
 
109
      CALL TBFGET (TID1,IDAT, FORM,LFIELD,DTYP,ISTAT)
 
110
      CALL TBFGET (TID1,IVAR, FORM,LFIELD,VTYP,ISTAT)
 
111
      CALL TBDGET (TID1,ISTORE,ISTAT)
 
112
      IF (ISTORE.NE.F_TRANS) THEN
 
113
        TEXT='Input table '//INAME1//' stored not transposed'
 
114
        CALL STETER(1,TEXT)
 
115
      ENDIF
 
116
      IF (TTYP.NE.D_R8_FORMAT.OR.DTYP.NE.D_R8_FORMAT.OR.
 
117
     $     VTYP.NE.D_R8_FORMAT) THEN
 
118
        CALL STETER(2,
 
119
     $    'Column(s) in 1st table must be in DOUBLE PRECISION')
 
120
      ENDIF
 
121
      CALL TBCMAP (TID1,ITIME,POBS1T,ISTAT)
 
122
      CALL TBCMAP (TID1,IDAT, POBS1D,ISTAT)
 
123
      CALL TBCMAP (TID1,IVAR, POBS1V,ISTAT)
 
124
C
 
125
C   Map input/output table
 
126
C
 
127
      CALL TBTOPN (INAME2,F_IO_MODE,TID2,ISTAT)
 
128
      CALL TBIGET (TID2,NCOL,NOBS2,ISOR,ICOL,NROW,ISTAT)
 
129
      CALL TBLSER (TID2,'TIME' ,ITIME,ISTAT)
 
130
      IF (ITIME.LT.0) THEN
 
131
        CALL STETER(8,'Column :TIME not found in 2nd table')
 
132
      ENDIF
 
133
      CALL TBLSER (TID2,'VALUE',IDAT ,ISTAT)
 
134
      IF (IDAT.LT.0) THEN
 
135
        CALL STETER(9,'Column :VALUE not found in 2nd table')
 
136
      ENDIF
 
137
      CALL TBLSER (TID2,'VAR'  ,IVAR, ISTAT)
 
138
      IF (IVAR.LT.0) THEN
 
139
        CALL STETER(10,'Column :VAR not found in 2nd table')
 
140
      ENDIF
 
141
      CALL TBFGET (TID2,ITIME,FORM,LFIELD,TTYP,ISTAT)
 
142
      CALL TBFGET (TID2,IDAT, FORM,LFIELD,DTYP,ISTAT)
 
143
      CALL TBFGET (TID2,IVAR, FORM,LFIELD,VTYP,ISTAT)
 
144
      CALL TBDGET (TID2,ISTORE,ISTAT)
 
145
      IF (ISTORE.NE.F_TRANS) THEN
 
146
        TEXT='Input table '//INAME2//' stored not transposed'
 
147
        CALL STETER(3,TEXT)
 
148
      ENDIF
 
149
      IF (TTYP.NE.D_R8_FORMAT.OR.DTYP.NE.D_R8_FORMAT.OR.
 
150
     $     VTYP.NE.D_R8_FORMAT) THEN
 
151
        CALL STETER(4,
 
152
     $    'Column(s) in 2nd table must be in DOUBLE PRECISION')
 
153
      ENDIF
 
154
      CALL TBCMAP (TID2,ITIME,POBS2T,ISTAT)
 
155
      CALL TBCMAP (TID2,IDAT, POBS2D,ISTAT)
 
156
      CALL TBCMAP (TID2,IVAR, POBS2V,ISTAT)
 
157
C
 
158
C   Create and map work tables  COBS,NU,CT
 
159
C
 
160
      ISIZE=(NOBS1+1)*(NOBS1+1)
 
161
      CALL STFCRE('ZZMIDAWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
 
162
     $   ISIZE,ICOBS,ISTAT)
 
163
      CALL STFMAP(ICOBS,F_X_MODE,1,ISIZE,ASIZE,PCOBS,ISTAT)
 
164
      ISIZE=NOBS1+1
 
165
      CALL STFCRE('ZZMIDBWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
 
166
     $   ISIZE,INU,  ISTAT)
 
167
      CALL STFMAP(INU,  F_X_MODE,1,ISIZE,ASIZE,PNU,ISTAT)
 
168
      ISIZE=NOBS1+1
 
169
      CALL STFCRE('ZZMIDCWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
 
170
     $   ISIZE,ICT,  ISTAT)
 
171
      CALL STFMAP(ICT,  F_X_MODE,1,ISIZE,ASIZE,PCT,ISTAT)
 
172
C
 
173
C   Interpolate using covariance function
 
174
C
 
175
      MODE=2
 
176
      IF (LLE(CFUNC,'LIN').AND.LGE(CFUNC,'LIN')) THEN
 
177
        CALL INTERPOL(
 
178
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
179
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
180
     $    ACFLIN,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
181
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
182
      ELSEIF (LLE(CFUNC,'POL').AND.LGE(CFUNC,'POL')) THEN
 
183
        CALL INTERPOL(
 
184
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
185
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
186
     $    ACFPOL,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
187
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
188
      ELSEIF (LLE(CFUNC,'POW').AND.LGE(CFUNC,'POW')) THEN
 
189
        CALL INTERPOL(
 
190
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
191
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
192
     $    ACFPOW,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
193
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
194
      ELSEIF (LLE(CFUNC,'EXP').AND.LGE(CFUNC,'EXP')) THEN
 
195
        CALL INTERPOL(
 
196
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
197
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
198
     $    ACFEXP,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
199
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
200
      ELSEIF (LLE(CFUNC,'LOG').AND.LGE(CFUNC,'LOG')) THEN
 
201
        CALL INTERPOL(
 
202
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
203
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
204
     $    ACFLOG,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
205
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
206
      ELSEIF (LLE(CFUNC,'IPO').AND.LGE(CFUNC,'IPO')) THEN
 
207
        CALL INTERPOL(
 
208
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
209
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
210
     $    ACFIPO,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
211
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
212
      ELSEIF (LLE(CFUNC,'UR0').AND.LGE(CFUNC,'UR0')) THEN
 
213
        CALL INTERPOL(
 
214
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
215
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
216
     $    TSADELUR0,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
217
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
218
      ELSEIF (LLE(CFUNC,'UR1').AND.LGE(CFUNC,'UR1')) THEN
 
219
        CALL INTERPOL(
 
220
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
221
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
222
     $    TSADELUR1,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
223
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
224
      ELSEIF (LLE(CFUNC,'UR2').AND.LGE(CFUNC,'UR2')) THEN
 
225
        CALL INTERPOL(
 
226
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
227
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
228
     $    TSADELUR2,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
229
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
230
      ELSEIF (LLE(CFUNC,'UR3').AND.LGE(CFUNC,'UR3')) THEN
 
231
        CALL INTERPOL(
 
232
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
233
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
234
     $    TSADELUR3,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
235
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
236
      ELSEIF (LLE(CFUNC,'UR4').AND.LGE(CFUNC,'UR4')) THEN
 
237
        CALL INTERPOL(
 
238
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
239
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
240
     $    TSADELUR4,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
241
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
242
      ELSEIF (LLE(CFUNC,'UR5').AND.LGE(CFUNC,'UR5')) THEN
 
243
        CALL INTERPOL(
 
244
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
245
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
246
     $    TSADELUR5,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
247
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
248
      ELSEIF (LLE(CFUNC,'UR6').AND.LGE(CFUNC,'UR6')) THEN
 
249
        CALL INTERPOL(
 
250
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
251
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
252
     $    TSADELUR6,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
253
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
254
      ELSEIF (LLE(CFUNC,'UR7').AND.LGE(CFUNC,'UR7')) THEN
 
255
        CALL INTERPOL(
 
256
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
257
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
258
     $    TSADELUR7,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
259
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
260
      ELSEIF (LLE(CFUNC,'UR8').AND.LGE(CFUNC,'UR8')) THEN
 
261
        CALL INTERPOL(
 
262
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
263
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
264
     $    TSADELUR8,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
265
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
266
      ELSEIF (LLE(CFUNC,'UR9').AND.LGE(CFUNC,'UR9')) THEN
 
267
        CALL INTERPOL(
 
268
     $    MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
 
269
     $    MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
 
270
     $    TSADELUR9,PARM,NOBS1,NOBS2,MODE,AVER,VAR,
 
271
     $    MADRID(PCOBS),MADRID(PNU),MADRID(PCT))
 
272
      ELSE
 
273
        CALL STETER(5,'Wrong name for ACF type')
 
274
      ENDIF
 
275
      NCOL=3
 
276
      CALL TBIPUT (TID2,NCOL,NOBS2,ISTAT)
 
277
C
 
278
C   Wind-up
 
279
C
 
280
      CALL DSCUPT(TID2,TID2,' ',ISTAT)
 
281
      CALL STSEPI
 
282
C
 
283
      END
 
284
C
 
285
C
 
286
C