← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/spec/src/integbin.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:
stdred/spec/src/integbin.for
 
1
C @(#)integbin.for      19.1 (ES0-DMD) 02/25/03 14:29:41
 
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 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding 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
      PROGRAM INGBIN
 
30
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
32
C                                         all rights reserved
 
33
C
 
34
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:13 - 3 DEC 1987
 
35
C
 
36
C.LANGUAGE: F77+ESOext
 
37
C
 
38
C.AUTHOR: D.PONZ
 
39
C
 
40
C.IDENTIFICATION
 
41
C  program    INTEGBIN              version 1.0          851028
 
42
C
 
43
C.Keywords
 
44
C  flux integration over fixed bins in 1-D images
 
45
C
 
46
C.PURPOSE
 
47
C   to be used for normalization of spectra or construction of 
 
48
C   response curves
 
49
C
 
50
C.ALGORITHM
 
51
C  get central positions of bins and their widths from table (in keyword
 
52
C  NFTAB), integrate flux in image (keyword NFIMA) over bins (where
 
53
C  necessary do a linear interpolation), write results to table
 
54
C
 
55
C.INPUT/OUTPUT
 
56
C the following keywords are used:
 
57
C NFIMA/C/1/60         name of input image
 
58
C NFTAB/C/1/60        name of table with central wavelengths, bin 
 
59
C                       widths, will also receive fluxes derived
 
60
C
 
61
C. modif MPeron 130291 add TBTCLO
 
62
C --------------------------------------------------------------------
 
63
C
 
64
      IMPLICIT NONE
 
65
      INTEGER MADRID
 
66
      CHARACTER*60 NFIMA,NFTAB
 
67
      CHARACTER    FORM*6,IDENT*72,CUNIT*72
 
68
      INTEGER   NPIX,ISTAT,IAV,ERROR(2),ICOL(3),IRCOL(2)
 
69
      INTEGER   NCOL,NROW,LL,I,NSC1
 
70
      INTEGER*8 PNTR
 
71
      INTEGER   NAXIS,KUN,KNUL
 
72
      INTEGER   IMNO,TID,NAC,NAR,DTYPE
 
73
      REAL      DATA(2),BINFLX
 
74
      DOUBLE PRECISION DEND, DSTEP, DSTART, LAMLO, LAMHI
 
75
      LOGICAL CURRY,NULL(3)
 
76
C
 
77
      INCLUDE  'MID_INCLUDE:TABLES.INC'
 
78
      COMMON /VMR/MADRID(1)
 
79
      INCLUDE  'MID_INCLUDE:TABLED.INC'
 
80
C
 
81
C ... set up MIDAS environment
 
82
C
 
83
      CALL STSPRO('INTEGBIN')
 
84
C
 
85
C ... get name of input frame, map it
 
86
C
 
87
      CALL STKRDC('NFIMAGE',1,1,60,IAV,NFIMA,KUN,KNUL,ISTAT)
 
88
      IF (ISTAT.NE.0) GO TO 20
 
89
      CALL STIGET(NFIMA,10,0,1,1,NAXIS,NPIX,DSTART,DSTEP,
 
90
     +IDENT,CUNIT,PNTR,IMNO,ISTAT)
 
91
      IF (ISTAT.NE.0) GO TO 20
 
92
      IF (NAXIS.NE.1) THEN
 
93
          CALL STTPUT('ERROR: input frame must be one-dimensional',
 
94
     +                ISTAT)
 
95
          ERROR(1) = 1
 
96
          GO TO 20
 
97
      END IF
 
98
      DEND    = DSTART + (NPIX-1)*DSTEP
 
99
C
 
100
C ... get name of table
 
101
C
 
102
      CALL STKRDC('NFTABLE',1,1,60,IAV,NFTAB,KUN,KNUL,ISTAT)
 
103
      IF (ISTAT.NE.0) GO TO 20
 
104
C
 
105
C ... intialize table
 
106
C
 
107
      CALL TBTOPN(NFTAB,2,TID,ISTAT)
 
108
      IF (ISTAT.NE.0) GO TO 20
 
109
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
 
110
      IF (ISTAT.NE.0) GO TO 20
 
111
C
 
112
C ... find column 'X_AXIS' (wavelength)
 
113
C
 
114
      CALL TBLSER(TID,'X_AXIS',ICOL(1),ISTAT)
 
115
      IF (ICOL(1).EQ.-1) THEN
 
116
          ISTAT  = 10
 
117
          CALL STTPUT('ERROR: column not found',ISTAT)
 
118
          GO TO 20
 
119
      END IF
 
120
      CALL TBFGET(TID,ICOL(1),FORM,LL,DTYPE,ISTAT)
 
121
      IF (DTYPE.EQ.D_C_FORMAT) THEN
 
122
          ISTAT  = 11
 
123
          CALL STTPUT('ERROR: wrong column format',ISTAT)
 
124
          GO TO 20
 
125
      END IF
 
126
C
 
127
C ... find column 'BIN_WIDTH'
 
128
C
 
129
      CALL TBLSER(TID,'BIN_WIDTH',ICOL(2),ISTAT)
 
130
      IF (ICOL(2).EQ.-1) THEN
 
131
          ISTAT  = 10
 
132
          CALL STTPUT('ERROR: column not found',ISTAT)
 
133
          GO TO 20
 
134
      END IF
 
135
      CALL TBFGET(TID,ICOL(2),FORM,LL,DTYPE,ISTAT)
 
136
      IF (DTYPE.EQ.D_C_FORMAT) THEN
 
137
          ISTAT  = 11
 
138
          CALL STTPUT('ERROR: wrong column format',ISTAT)
 
139
          GO TO 20
 
140
      END IF
 
141
C
 
142
C ... find column 'Y_AXIS' (flux)
 
143
C
 
144
      CALL TBLSER(TID,'Y_AXIS',ICOL(3),ISTAT)
 
145
      IF (ICOL(3).EQ.-1) THEN
 
146
          ISTAT  = 10
 
147
          CALL STTPUT('ERROR: column not found',ISTAT)
 
148
          GO TO 20
 
149
      END IF
 
150
      CALL TBFGET(TID,ICOL(3),FORM,LL,DTYPE,ISTAT)
 
151
      IF (DTYPE.EQ.D_C_FORMAT) THEN
 
152
          ISTAT  = 11
 
153
          CALL STTPUT('ERROR: wrong column format',ISTAT)
 
154
          GO TO 20
 
155
      END IF
 
156
      IRCOL(1) = ICOL(1)
 
157
      IRCOL(2) = ICOL(2)
 
158
C
 
159
C ... read table
 
160
C
 
161
      DO 10 I = 1,NROW
 
162
          CALL TBRRDR(TID,I,2,IRCOL,DATA,NULL,ISTAT)
 
163
          CURRY = ( .NOT. NULL(1)) .AND. ( .NOT. NULL(2))
 
164
          IF (CURRY) THEN  !   lower limit in lambda of bin
 
165
              LAMLO  = DATA(1) - DATA(2)/2.  !   upper limit in lambda of bin
 
166
              LAMHI  = DATA(1) + DATA(2)/2.
 
167
              IF ((LAMLO.GE.DSTART) .AND. (LAMHI.LE.DEND)) THEN
 
168
C
 
169
C ... let subroutine DOINTG do the job
 
170
C
 
171
                CALL DOINTG(MADRID(PNTR),DSTART,DSTEP,NPIX,LAMLO,LAMHI,
 
172
     +                         BINFLX)
 
173
                CALL TBEWRR(TID,I,ICOL(3),BINFLX,ISTAT)
 
174
              END IF
 
175
          END IF
 
176
   10 CONTINUE
 
177
C
 
178
C the end
 
179
C
 
180
      CALL TBTCLO(TID,ISTAT)
 
181
      CALL STSEPI
 
182
      STOP 
 
183
C
 
184
C ... direct trouble makers here:
 
185
C
 
186
   20 ERROR(1) = ISTAT
 
187
      CALL STKWRI('PROGSTAT',ERROR,1,2,KUN,ISTAT)
 
188
      CALL STSEPI
 
189
      STOP 
 
190
C
 
191
      END
 
192
 
 
193
      SUBROUTINE DOINTG(FLUX,DSTART,DSTEP,NPIX,LAMLO,LAMHI,BINFLX)
 
194
C
 
195
C
 
196
      IMPLICIT NONE
 
197
      INTEGER NPIX,IPIXHI,IPIXLO,I
 
198
      REAL FLUX(1),PIXLO,PIXHI,FPIXHI,FPIXLO,BINFLX
 
199
      DOUBLE PRECISION DSTART, DSTEP, LAMLO, LAMHI
 
200
C
 
201
C ... all operations below assume that the flux is evenly distributed
 
202
C ... ver the pixel and that its coordinates refer to the pixel's 
 
203
C ... center !
 
204
C  !   first pixel hit (only partly)
 
205
      PIXLO  = 0.5 + (LAMLO-DSTART)/DSTEP  !   last pixel hit (also only partly)
 
206
      PIXHI  = 0.5 + (LAMHI-DSTART)/DSTEP  !   nmubers of first ...
 
207
      IPIXLO = INT(PIXLO)  !   ... and last pixel
 
208
      IPIXHI = INT(PIXHI)  !   fraction of first and ...
 
209
      FPIXLO = IPIXLO + 1 - PIXLO  !   ... last pixel to be considered
 
210
      FPIXHI = PIXHI - IPIXHI
 
211
C
 
212
C ... do a linear interpolation across the first and last pixel,
 
213
C ... add up the pixels in between, return mean flux per pixel
 
214
C
 
215
      BINFLX = FPIXLO*FLUX(IPIXLO) + FPIXHI*FLUX(IPIXHI)
 
216
      DO 10 I = IPIXLO + 1,IPIXHI - 1
 
217
          BINFLX = BINFLX + FLUX(I)
 
218
   10 CONTINUE
 
219
      BINFLX = BINFLX/ (IPIXHI-IPIXLO-1+FPIXLO+FPIXHI)
 
220
C
 
221
      RETURN
 
222
      END