← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necrebi.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/echelle/src/necrebi.for
 
1
C @(#)necrebi.for       19.1 (ESO-IPG) 02/25/03 14:20: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
      PROGRAM ECHREB
 
30
C
 
31
C+
 
32
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
33
C                                         all rights reserved
 
34
C
 
35
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  10:37 - 11 JAN 1988
 
36
C
 
37
C.LANGUAGE: F77+ESOext
 
38
C
 
39
C.AUTHOR: D.PONZ
 
40
C
 
41
C
 
42
C.IDENTIFICATION
 
43
C
 
44
C  program  ECHREBI1
 
45
C
 
46
C.MODIFICATIONS
 
47
C
 
48
C  910128   P.Ballester   Define and set variables
 
49
C
 
50
C.KEYWORDS
 
51
C
 
52
C  ECHELLE, CASPEC, REBIN
 
53
C
 
54
C.PURPOSE
 
55
C
 
56
C  OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM
 
57
C  OF TYPE ECHELLE WITH THE DISPERSION COEFFS DEFINED AS DESCRIPTORS
 
58
C
 
59
C.ALGORITHM
 
60
C
 
61
C  EACH LINE IN THE INPUT FRM IS REBINNED IN LINEAR STEP IN WAVELENGTH
 
62
C  ACCORDING TO THE DISPERSION COEFFS.
 
63
C
 
64
C.INPUT/OUTPUT
 
65
C
 
66
C  REBIN/ECHELLE INPUT OUTPUT STEP [METHOD] [TABLE]
 
67
C  OR
 
68
C  REBIN/ECHELLE INPUT OUTPUT REFER_IMAGE [METHOD] [TABLE]
 
69
C
 
70
C METHOD CAN BE EITHER "NONLINEAR" - NON LINEAR REBINNING USING
 
71
C                                     THE COEFFS IN THE TABLE
 
72
C                      "LINEAR"    - LINEAR REBINNING
 
73
C
 
74
C.MODIFICATIONS
 
75
C [2.2]  DEFINE THE SAMPLING DOMAIN VIA A REFER IMAGE
 
76
C 
 
77
C 010201                last modif
 
78
C 
 
79
C 
 
80
C-
 
81
      IMPLICIT NONE
 
82
C
 
83
      INTEGER      MADRID, ID, IR
 
84
      INTEGER      I,I1,IMNO
 
85
      INTEGER*8    IPNTRA
 
86
      INTEGER*8    IPNTRB
 
87
      INTEGER      ISTAT,NAXISA,NAXISB
 
88
      INTEGER      PIXBN, KUN,KNUL, NNNN
 
89
      INTEGER      NPIXA(3),NPIXB(3),IORD(2),IORD2(2),TID
 
90
      INTEGER      IORDER(100),NCOE(2),IA1(100),NPTOT(100)
 
91
      DOUBLE PRECISION WSTART(100),STEPA(3),STARTA(3)
 
92
      DOUBLE PRECISION STEPB(3),STARTB(3)
 
93
      DOUBLE PRECISION WSTEP, A1(7,100)
 
94
      REAL         CUT(4),XMIN,XMAX
 
95
      CHARACTER*60 OUTFRM,INFRM,REFER,TABLE
 
96
      CHARACTER    METHOD*1, LINE*80, DSCNAM*8
 
97
      CHARACTER    IDENTA*72,CUNITA*64,CUNITB*64
 
98
      LOGICAL      IREF
 
99
      INCLUDE      'MID_INCLUDE:ST_DEF.INC'
 
100
      COMMON /VMR/MADRID(1)
 
101
      INCLUDE      'MID_INCLUDE:ST_DAT.INC'
 
102
      DATA CUNITB/'Flux            Relative wavel.'/
 
103
      DATA NPIXA /3*1/, STARTA /3*0.0/, STEPA /3*1.0/
 
104
      DATA NPIXB /3*1/, STARTB /3*0.0/, STEPB /3*1.0/
 
105
C
 
106
C ... INITIALIZE SYSTEM
 
107
C
 
108
      CALL STSPRO('ECHREB')
 
109
      CALL STKRDC('IN_A',1,1,60,I,INFRM,KUN,KNUL,ISTAT)
 
110
      CALL STKRDC('OUT_A',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT)
 
111
      CALL CLNFRA(INFRM,INFRM,0)
 
112
      CALL CLNFRA(OUTFRM,OUTFRM,0)
 
113
      CALL STKRDC('P3',1,1,60,I,REFER,KUN,KNUL,ISTAT)
 
114
      CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT)
 
115
      I1     = INDEX('.+1234567890',REFER(1:1))
 
116
      IF (I1.NE.0) THEN
 
117
          CALL GENCNV(REFER,4,1,NNNN,XMIN,WSTEP,ISTAT)
 
118
          IF (ISTAT.LT.1) CALL STETER(11,'Wrong sampling step')
 
119
C 
 
120
          WRITE(LINE,9010) WSTEP
 
121
          CALL STTPUT(LINE,ISTAT)
 
122
          IREF   = .FALSE.
 
123
      ELSE
 
124
          IREF   = .TRUE.
 
125
          CALL CLNFRA(REFER,REFER,0)
 
126
      END IF
 
127
      IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN
 
128
          CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,ISTAT)
 
129
          CALL CLNTAB(TABLE,TABLE,0)
 
130
          CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
 
131
          CALL STDRDI(TID,'ORDER',1,2,I,IORD,KUN,KNUL,ISTAT)
 
132
          CALL STDRDI(TID,'COEFS',1,2,I,NCOE,KUN,KNUL,ISTAT)
 
133
          CALL STDRDD(TID,'COEFD',1,7*NCOE(2),I,A1,
 
134
     .                KUN,KNUL,ISTAT)
 
135
          CALL STDRDI(TID,'COEFI',1,NCOE(2),I,IA1,
 
136
     .                KUN,KNUL,ISTAT)
 
137
          CALL TBTCLO(TID,ISTAT)
 
138
      END IF
 
139
C
 
140
C ... MAP INPUT IMAGE
 
141
C
 
142
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
143
     .            3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
 
144
     .            IPNTRA,IMNO,ISTAT)
 
145
 
 
146
      IF (NAXISA.EQ.1 .OR. NPIXA(2).EQ.1) THEN
 
147
          NPIXA(2) = 1
 
148
          NAXISA = 1
 
149
      END IF
 
150
 
 
151
      IF (NPIXA(2).GT.100)  THEN
 
152
          CALL STETER(9,
 
153
     +   'Module ECHRE31: Too many orders (limit 100). Buffer overflow')
 
154
      ENDIF
 
155
 
 
156
      IF (STEPA(2).NE.1.)  THEN
 
157
          CALL STETER(9,'Module ECHREBI: Wrong step(2) .ne. 1 ')
 
158
      ENDIF
 
159
 
 
160
      IORD2(1) = STARTA(2)                  ! Starting relative order number
 
161
      IORD2(2) = STARTA(2) + NPIXA(2) - 1   ! Final relative order number
 
162
C
 
163
C ... MAP OUTPUT FRM
 
164
C
 
165
      IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN
 
166
          NAXISB = 2
 
167
          IF (IREF) THEN
 
168
              CALL STFOPN(REFER,D_R4_FORMAT,0,F_IMA_TYPE,
 
169
     .                    IR,ISTAT)
 
170
              CALL STDRDI(IR,'NPIX',1,2,I,NPIXB,KUN,KNUL,ISTAT)
 
171
              CALL STDRDD(IR,'START',1,2,I,STARTB,KUN,KNUL,ISTAT)
 
172
              CALL STDRDD(IR,'STEP',1,2,I,STEPB,KUN,KNUL,ISTAT)
 
173
              WSTEP  = STEPB(1)
 
174
              CALL STDRDI(IR,'NORDER',1,NPIXB(2),I,IORDER,
 
175
     .                    KUN,KNUL,ISTAT)
 
176
              CALL STDRDI(IR,'NPTOT',1,NPIXB(2),I,NPTOT,
 
177
     .                    KUN,KNUL,ISTAT)
 
178
              CALL STDRDD(IR,'WSTART',1,NPIXB(2),I,WSTART,
 
179
     .                    KUN,KNUL,ISTAT)
 
180
              DO 10 I = 1,NPIXB(2)
 
181
                  WSTART(I) = WSTART(I)/1000.
 
182
   10         CONTINUE
 
183
          ELSE
 
184
              CALL WRANGE(NCOE(2),A1,IA1,WSTEP,IORD,NPIXA(1),WSTART,
 
185
     .                    NPIXB(1),NPTOT,IORDER,STARTA,STEPA,IORD2)
 
186
              NPIXB(2)  = NPIXA(2)
 
187
              STARTB(1) = 0.D0
 
188
              STARTB(2) = STARTA(2)
 
189
              STEPB(1)  = WSTEP
 
190
              STEPB(2)  = STEPA(2)
 
191
          END IF
 
192
 
 
193
          CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
194
     .                NAXISB,NPIXB,STARTB,STEPB,
 
195
     .                IDENTA,CUNITB,IPNTRB,ID,ISTAT)
 
196
C
 
197
C ... NON LINEAR REBINNING
 
198
C
 
199
          CALL ECHRE3(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA,
 
200
     .                NCOE(2),A1,IA1,IORD,MADRID(IPNTRB),NPIXB(1),
 
201
     .                NPIXB(2),WSTART,WSTEP,IORD2)
 
202
C
 
203
C ... WRITE PROCESS DESCRIPTORS
 
204
C
 
205
          NNNN = NPIXB(1)*NPIXB(2)
 
206
          CALL MNMX(MADRID(IPNTRB),NNNN,XMIN,XMAX)
 
207
          CUT(1) = XMIN
 
208
          CUT(2) = XMAX
 
209
          CUT(3) = XMIN
 
210
          CUT(4) = XMAX
 
211
          DO 20 I = 1,NPIXB(2)
 
212
              WSTART(I) = 1000.*WSTART(I)
 
213
   20     CONTINUE
 
214
          CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT)
 
215
          CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT)
 
216
          CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT)
 
217
          CALL STDWRI(ID,'NORDER',IORDER,1,NPIXB(2),KUN,ISTAT)
 
218
          CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT)
 
219
      ELSE
 
220
C ... LINEAR REBINNING
 
221
          IF (NPIXA(2).GT.1) THEN
 
222
CCCCCCCCCCCCC              CALL STDFND(IMNO,'WSTART',I1,I2,I3,ISTAT)
 
223
C              IF (TYPE(1:1).EQ.' ') THEN
 
224
C                  CALL STTPUT(' Wrong input image ... ',ISTAT)
 
225
C                  GO TO 40
 
226
CCCCCCCCCCCCC              END IF
 
227
              CALL STDRDD(IMNO,'WSTART',1,NPIXA(2),I,WSTART,
 
228
     .                    KUN,KNUL,ISTAT)
 
229
              CALL STDRDI(IMNO,'NPTOT',1,NPIXA(2),I,NPTOT,
 
230
     .                    KUN,KNUL,ISTAT)
 
231
          END IF
 
232
          CALL STDRDR(IMNO,'LHCUTS',1,4,I,CUT,KUN,KNUL,ISTAT)
 
233
          PIXBN     = STEPA(1)/WSTEP
 
234
          NAXISB    = NAXISA
 
235
          NPIXB(1)  = NPIXA(1)*PIXBN
 
236
          NPIXB(2)  = NPIXA(2)
 
237
          STARTB(1) = STARTA(1)
 
238
          STARTB(2) = STARTA(2)
 
239
          STEPB(1)  = WSTEP
 
240
          STEPB(2)  = STEPA(2)
 
241
          IF (NPIXB(2).GT.1) THEN
 
242
              DO 30 I = 1,NPIXB(2)
 
243
                  NPTOT(I) = NPTOT(I)*PIXBN
 
244
   30         CONTINUE
 
245
          ELSE
 
246
              WSTART(1) = STARTB(1)
 
247
              NPTOT(1) = NPIXB(1)
 
248
          END IF
 
249
          CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
250
     .                NAXISB,NPIXB,STARTB,STEPB,IDENTA,
 
251
     .                CUNITB,IPNTRB,ID,ISTAT)
 
252
          CALL ECHR01(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA,
 
253
     .                NPTOT,WSTART,
 
254
     .                MADRID(IPNTRB),NPIXB(1),NPIXB(2),WSTEP)
 
255
 
 
256
          CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT)
 
257
          IF (NPIXB(2).GT.1) THEN
 
258
            CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT)
 
259
            CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT)
 
260
          END IF
 
261
          CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT)
 
262
      END IF
 
263
      CALL STSEPI
 
264
C10000 FORMAT(G<I1>.0)
 
265
 9000 FORMAT (G6.0)
 
266
 9010 FORMAT(' Sampling step = ',F10.4)
 
267
      END
 
268
 
 
269
      SUBROUTINE ECHRE3(X,NPIXA1,NPIXA2,STARTA,STEPA,NO,A,IA,IORD,
 
270
     +                  Y,NPIXB1,NPIXB2,YSTR,YSTP,IORD2)
 
271
C
 
272
C REBIN THE ORDERS
 
273
C
 
274
      IMPLICIT NONE
 
275
C
 
276
      INTEGER NO, NPIXB1, NPIXB2
 
277
      INTEGER NPIXA1,NPIXA2,IORD(2),IA(NO),IORD2(2)
 
278
      REAL X(NPIXA1,NPIXA2)
 
279
      REAL Y(NPIXB1,NPIXB2)
 
280
      DOUBLE PRECISION YSTR(1), YSTP, STARTA(2), STEPA(2)
 
281
      DOUBLE PRECISION A(7,NO),B1(3)
 
282
C
 
283
      INTEGER KMIN, KMAX, I0, I2, I1, K1, II, MIN, MAX, I3
 
284
      DOUBLE PRECISION WLSTP, WLST2
 
285
C
 
286
C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS)
 
287
C
 
288
      KMIN   = MIN(IORD(1),IORD(2))
 
289
      KMAX   = MAX(IORD(1),IORD(2))
 
290
      IF (IORD(1).LT.IORD(2)) THEN
 
291
          I0     = NO
 
292
          I2     = -1
 
293
      ELSE
 
294
          I0     = 1  
 
295
          I2     = 1
 
296
      END IF
 
297
      I1     = I0
 
298
      WLSTP  = YSTP/1000.
 
299
C
 
300
C ... ITERATION ON ORDERS
 
301
C
 
302
      DO 20 K1 = IORD2(2), IORD2(1), -1
 
303
          I3 = I1 + IORD2(1) - 1
 
304
          IF (IA(I3).LT.0) THEN
 
305
              B1(1)  = A(1,I3)
 
306
              B1(2)  = -A(2,I3)
 
307
              B1(3)  = A(3,I3)
 
308
          END IF
 
309
          WLST2  = YSTR(I1)
 
310
          IF (IA(I3).LT.0) THEN
 
311
              CALL LREBIN(X(1,I1),NPIXA1,STARTA,STEPA,B1,3,Y(1,I1),
 
312
     +                    NPIXB1,WLST2,WLSTP)
 
313
C
 
314
C ... SIGN CHANGE
 
315
              DO 10 II = 1,NPIXB1
 
316
                  Y(II,I1) = -Y(II,I1)
 
317
   10         CONTINUE
 
318
          ELSE
 
319
              CALL LREBI1(X(1,I1),NPIXA1,STARTA,STEPA,A(1,I3),IA(I3),
 
320
     +                     Y(1,I1),NPIXB1,WLST2,WLSTP)
 
321
          END IF
 
322
          I1     = I1 + I2
 
323
          IF (I1.GT.NO) RETURN
 
324
   20 CONTINUE
 
325
      RETURN
 
326
      END
 
327
 
 
328
      DOUBLE PRECISION FUNCTION POLEV(IA,A,XX)
 
329
C
 
330
C
 
331
C
 
332
      IMPLICIT NONE
 
333
C
 
334
      INTEGER IA
 
335
      DOUBLE PRECISION A(IA),XX,YY
 
336
 
 
337
      YY     = A(1) + XX* (A(2)+XX* (A(3)+XX* (A(4)+XX* (A(5)+XX* (A(6)+
 
338
     +         XX*A(7))))))
 
339
      POLEV  = YY
 
340
      RETURN
 
341
      END
 
342
 
 
343
      SUBROUTINE MNMX(X,N,XMIN,XMAX)
 
344
C
 
345
C FIND MIN MAX VALS
 
346
C
 
347
      IMPLICIT NONE
 
348
      INTEGER  N
 
349
      REAL     X(N), XMIN, XMAX
 
350
C
 
351
      INTEGER I
 
352
C
 
353
      XMIN   = X(1)
 
354
      XMAX   = XMIN
 
355
      DO 10 I = 2,N
 
356
          XMIN   = AMIN1(XMIN,X(I))
 
357
          XMAX   = AMAX1(XMAX,X(I))
 
358
   10 CONTINUE
 
359
      RETURN
 
360
      END
 
361
 
 
362
      SUBROUTINE WRANGE(NCOE,A,IA,WSTEP,IORD,NPIXA,WSTART,NPIXB,NPTOT,
 
363
     +                  IORDER,STARTA,STEPA,IORD2)
 
364
C
 
365
C GET STARTING POSITIONS FOR EACH ORDER AND NO. OF PIXELS PER ORDER
 
366
C
 
367
C
 
368
C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS)
 
369
C
 
370
      IMPLICIT NONE
 
371
C
 
372
      INTEGER NCOE, NPIXA, NPIXB
 
373
      DOUBLE PRECISION A(7,NCOE),POLEV,B(3), STARTA(1), STEPA(1)
 
374
      DOUBLE PRECISION WSTART(1), WSTEP, WLST1, WLEN1
 
375
      INTEGER IA(NCOE),IORD(2),NPTOT(1),IORDER(1), IORD2(2)
 
376
C
 
377
      INTEGER KMIN, KMAX, I0, I1, I2, I, ISTEP, K, N1, I3
 
378
c      INTEGER*8 N1
 
379
      DOUBLE PRECISION WS1, WW1, IX_IR8
 
380
C
 
381
      KMIN   = MIN(IORD(1),IORD(2))
 
382
      KMAX   = MAX(IORD(1),IORD(2))
 
383
      IF (IORD(1).LT.IORD(2)) THEN
 
384
          I0     = NCOE 
 
385
          I2     = -1
 
386
      ELSE
 
387
          I0     = 1 
 
388
          I2     = 1
 
389
      END IF
 
390
      I1     = I0
 
391
      NPIXB  = 0
 
392
      WS1    = WSTEP/1000.
 
393
C
 
394
C ... get order number
 
395
C
 
396
      ISTEP  = 1
 
397
      IF (IORD(1).GT.IORD(2)) ISTEP  = -1 ! Default case in present versions
 
398
      IORDER(1) = IORD(1) - IORD2(1) + 1
 
399
      DO 10 I = 2, IORD2(2) - IORD2(1) + 1
 
400
          IORDER(I) = IORDER(I-1) + ISTEP
 
401
   10 CONTINUE
 
402
C
 
403
C ... iteration on orders to get starting wavelength and no. of pixels
 
404
C
 
405
      DO 20 K = IORD2(2), IORD2(1), -1
 
406
          I3 = I1 + IORD2(1) - 1
 
407
          IF (IA(I3).LT.0) THEN
 
408
              B(1)   = A(1,I3)
 
409
              B(2)   = -A(2,I3)
 
410
              B(3)   = A(3,I3)
 
411
              CALL ECHORD(B,3,NPIXA,WLST1,WLEN1,STARTA,STEPA)
 
412
          ELSE
 
413
              WW1    = POLEV(IA(I3),A(1,I3),IX_IR8(1,STARTA,STEPA))
 
414
              WLST1  = WW1
 
415
              WW1    = POLEV(IA(I3),A(1,I3),IX_IR8(NPIXA,STARTA,STEPA))
 
416
              WLEN1  = WW1
 
417
          END IF
 
418
          WLST1  = WLST1*10.
 
419
          WLEN1  = WLEN1*10.
 
420
          N1     = (WLEN1-WLST1)/WS1
 
421
          IF (N1.GT.65534) THEN !added swolf@eso.org
 
422
             N1 = 0             !if WLEN1 and WLST1 are NaN N1 gets wrong!
 
423
          ENDIF                 !--> NPIXB gets wrong, too!
 
424
          NPTOT(I1) = N1
 
425
          NPIXB  = MAX(NPIXB,N1)
 
426
          IF (I1.EQ.I0) THEN
 
427
              WSTART(I1) = WS1*NINT(WLST1/WS1) !before: WLST1 (swolf@eso.org)
 
428
          ELSE
 
429
              WSTART(I1) = WSTART(I0) + WS1*NINT((WLST1-WSTART(I0))/WS1)
 
430
          END IF
 
431
          I1     = I1 + I2
 
432
          IF (I1.GT.NCOE) GO TO 30
 
433
   20 CONTINUE
 
434
   30 NPIXB  = NPIXB + 3
 
435
      RETURN
 
436
      END
 
437
 
 
438
      SUBROUTINE ECHR01(A,NPIXA1,NPIXA2,STARTA,STEPA,NPTOT,WSTART,
 
439
     +                  B,NPIXB1,NPIXB2,WSTEP)
 
440
C
 
441
C REBIN ECHELLE ORDERS
 
442
C
 
443
      IMPLICIT NONE
 
444
C
 
445
      INTEGER NPIXA1,NPIXA2,NPIXB1,NPIXB2
 
446
      DOUBLE PRECISION  STEPA(2),STARTA(2)
 
447
      REAL A(NPIXA1,NPIXA2),B(NPIXB1,NPIXB2)
 
448
      DOUBLE PRECISION  WSTART(NPIXA2), WSTEP
 
449
      INTEGER NPTOT(NPIXA2)
 
450
C
 
451
      INTEGER NORD, IORD, I, NPREB
 
452
      DOUBLE PRECISION XSTART, XP, ADEX, PIXBN
 
453
      REAL   VBIN
 
454
C
 
455
      PIXBN  = WSTEP/STEPA(1)
 
456
      NORD   = NPIXA2
 
457
      DO 30 IORD = 1,NORD
 
458
          XSTART = WSTART(IORD)
 
459
          NPREB  = NPTOT(IORD)
 
460
          DO 10 I = 1,NPREB
 
461
              XP     = XSTART + (I-1)*WSTEP
 
462
              ADEX   = (XP-XSTART)/STEPA(1) + 1.
 
463
              B(I,IORD) = VBIN(A(1,IORD),PIXBN,ADEX)
 
464
   10     CONTINUE
 
465
          DO 20 I = NPREB + 1,NPIXB1
 
466
              B(I,IORD) = 0.
 
467
   20     CONTINUE
 
468
   30 CONTINUE
 
469
      RETURN
 
470
      END
 
471
 
 
472
      REAL FUNCTION VBIN(ARRAY,PIXPBN,ADEX)
 
473
C
 
474
C FUNCTION TO RETURN A VALUE FOR A REBINNED BIN
 
475
C
 
476
C INPUT ARGUMENTS
 
477
C ARRAY REAL*4 ARRAY WITH REAL ORIGINAL PIXELS
 
478
C PIXPBN REAL*4 PIXELS PER BIN, NUMBER OF ORIGINAL
 
479
C   PIXELS IN THE NEW BIN, WITH FRACTIONS
 
480
C ADEX REAL*4 CENTER INDEX IN THE ARRAY FOR THE NEW BIN,
 
481
C   WITH FRACTIONS
 
482
C
 
483
      IMPLICIT NONE
 
484
      INTEGER          IBOTDX,ITOPDX,IENLP,I,ISTLP
 
485
      REAL             HBIN,BOTDX,TOPDX,SUM,BOTFRA,TOPFRA
 
486
      REAL             ARRAY(1)
 
487
      DOUBLE PRECISION PIXPBN, ADEX
 
488
C  !   HALF BIN
 
489
      HBIN   = PIXPBN*0.5  !   BOTTOM POINTER, FRACTION
 
490
      BOTDX  = ADEX - HBIN  !   TOP POINTER, FRACTION
 
491
      TOPDX  = ADEX + HBIN  !   BOTTOM POINTER, WITHOUT FRACTION
 
492
      IBOTDX = BOTDX  !   TOP POINTER, WITHOUT FRACTION
 
493
      ITOPDX = TOPDX
 
494
C
 
495
      IF (IBOTDX.EQ.ITOPDX) THEN
 
496
C
 
497
C   BIN SMALLER THAN ORIGINAL PIXELS
 
498
C   AND NEW BIN WITHIN AN ORIGINAL PIXEL
 
499
C
 
500
          SUM    = PIXPBN*ARRAY(IBOTDX+1)
 
501
      ELSE
 
502
C
 
503
C   BIN LARGER THAN ORIGINAL PIXELS, OR
 
504
C   SMALL BIN OVER TWO PIXELS
 
505
C
 
506
          BOTFRA = 1. - (BOTDX-FLOAT(IBOTDX))
 
507
          TOPFRA = TOPDX - FLOAT(ITOPDX)  !   START LOOP INDEX
 
508
          ISTLP  = IBOTDX + 2
 
509
          IENLP  = ITOPDX
 
510
          IBOTDX = IBOTDX + 1
 
511
          ITOPDX = ITOPDX + 1
 
512
C
 
513
          SUM    = 0.
 
514
          IF (IENLP.GE.ISTLP) THEN
 
515
              DO 10 I = ISTLP,IENLP
 
516
                  SUM    = SUM + ARRAY(I)
 
517
   10         CONTINUE
 
518
          END IF  !   LOWER FRACTION
 
519
          SUM    = SUM + BOTFRA*ARRAY(IBOTDX)  !   UPPER FRACTION
 
520
          SUM    = SUM + TOPFRA*ARRAY(ITOPDX)
 
521
      END IF
 
522
C
 
523
C AREA CORRECTION
 
524
C
 
525
      VBIN   = SUM/PIXPBN
 
526
      RETURN
 
527
      END
 
528
 
 
529
      SUBROUTINE ECHORD(AC,NC,NPIX,W1,W2,STARTA,STEPA)
 
530
C
 
531
C
 
532
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
533
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
534
C                                         all rights reserved
 
535
C
 
536
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  22:34 - 3 DEC 1987
 
537
C
 
538
C.LANGUAGE: F77+ESOext
 
539
C
 
540
C.AUTHOR: D.PONZ
 
541
C
 
542
C.KEYWORDS:
 
543
C  SPECTROSCOPY, CASPEC, ORDER DEFINITION
 
544
C
 
545
C.PURPOSE:
 
546
C  GET WAVELENGTH RANGE OF A GIVEN ORDER
 
547
C
 
548
C.ALGORITHM:
 
549
C
 
550
C.INPUT/OUTPUT
 
551
C  INPUT ARGUMENTS
 
552
C
 
553
C  AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION
 
554
C  NC  REAL*4 NO. OF COEFFICIENTS
 
555
C  NPIX  INTG*4 NO OF SAMPLES
 
556
C
 
557
C  OUTPUT ARGUMENTS
 
558
C
 
559
C  W1  REAL*4  LOWER LIMIT
 
560
C  W2  REAL*4  UPPER LIMIT
 
561
C
 
562
C-------------------------------------------------------------
 
563
C
 
564
C
 
565
      IMPLICIT NONE
 
566
      INTEGER NPIX, NC
 
567
      DOUBLE PRECISION AC(NC), STARTA(1), STEPA(1)
 
568
      DOUBLE PRECISION W1, W2, IX_IR8
 
569
C
 
570
      DOUBLE PRECISION A0,A1,A2,F1,FN
 
571
C
 
572
C
 
573
C ... ASSIGN COEFFS
 
574
C
 
575
C      F1     = 1.D0
 
576
C      FN     = NPIX
 
577
 
 
578
      F1     = IX_IR8(1,STARTA(1),STEPA(1))
 
579
      FN     = IX_IR8(NPIX,STARTA(1),STEPA(1))
 
580
 
 
581
      A0     = AC(1)
 
582
      A1     = AC(2)
 
583
      IF (NC.GE.3) THEN
 
584
          A2     = AC(3)
 
585
          W1     = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-F1)))/ (2.D0*A2)
 
586
          W2     = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-FN)))/ (2.D0*A2)
 
587
 
 
588
      ELSE
 
589
          W1     = (F1-A0)/A1
 
590
          W2     = (FN-A0)/A1
 
591
      END IF
 
592
      RETURN
 
593
      END
 
594
 
 
595
 
 
596
      SUBROUTINE LREBI1(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP)
 
597
C
 
598
C
 
599
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
600
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
601
C                                         all rights reserved
 
602
C
 
603
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  0:22 - 4 DEC 1987
 
604
C
 
605
C.LANGUAGE: F77+ESOext
 
606
C
 
607
C.AUTHOR: D.PONZ
 
608
C
 
609
C.IDENTIFICATION
 
610
C  SUBROUTINE  LREBIN1 VERSION 1   22 MAR 1983
 
611
C  J.D.PONZ
 
612
C
 
613
C.KEYWORDS:
 
614
C  SPECTROSCOPY, CASPEC, REBINNING
 
615
C
 
616
C.PURPOSE:
 
617
C  DATA REBINNING
 
618
C  CONVERSION FROM SAMPLES INTO WAVELENGTHS
 
619
C
 
620
C.ALGORITHM:
 
621
C  LINEAR INTERPOLATION. THE TRANSFORMATION SAMPLE NUMBER INTO
 
622
C  WAVELENGTH IS DEFINED BY
 
623
C  W = A0 + A1*S + A2*S**2 + ...
 
624
C
 
625
C.INPUT/OUTPUT
 
626
C  INPUT ARGUMENTS
 
627
C
 
628
C  X(NX) REAL*4 INPUT ARRAY
 
629
C  NX  INTG*4 NO OF SAMPLES
 
630
C  XSTR  REAL*4 STARTING VALUE OF INPUT
 
631
C  XSTP  REAL*4 INPUT STEP
 
632
C  AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION
 
633
C  NC  REAL*4 NO. OF COEFFICIENTS
 
634
C  NY  INTG*4 NO. OF OUTPUT VALUES
 
635
C  YSTR  REAL*4 STARTING VALUE OF THE INPUT
 
636
C  YSTP  REAL*4 STEP OF THE OUTPUT
 
637
C
 
638
C  OUTPUT ARGUMENTS
 
639
C
 
640
C  Y(NY) REAL*4 OUTPUT ARRAY
 
641
C
 
642
C-------------------------------------------------------------
 
643
C
 
644
C
 
645
      IMPLICIT NONE
 
646
      INTEGER  NX, NY, NC
 
647
      REAL X(NX),Y(NY)
 
648
      DOUBLE PRECISION YSTR, YSTP, XSTR, XSTP
 
649
C
 
650
      DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,A6,FN,FL,WPOLY
 
651
      DOUBLE PRECISION FACTOR,DPOLY,CNORM, IX_R8R8, XN, XNP1
 
652
      DOUBLE PRECISION PN,PNP1,TM,TMP1,C,COV,VAL, COVX
 
653
      DOUBLE PRECISION B2, B3, B4, B5, B6, DELT, STARTL
 
654
      INTEGER          M1, MST, N, K, MPT
 
655
C
 
656
      WPOLY(FL) = A0 + FL* (A1+FL* (A2+FL* (A3+FL* (A4+FL* (A5+FL*
 
657
     +            A6)))))
 
658
      DPOLY(FL) = A1 + FL* (B2+FL* (B3+FL* (B4+FL* (B5+FL*B6))))
 
659
C
 
660
C ... ASSIGN COEFFS
 
661
C
 
662
      A0     = AC(1)
 
663
      A1     = AC(2)
 
664
      A2     = AC(3)
 
665
      A3     = AC(4)
 
666
      A4     = AC(5)
 
667
      A5     = AC(6)
 
668
      A6     = AC(7)
 
669
      B2     = 2.D0*A2
 
670
      B3     = 3.D0*A3
 
671
      B4     = 4.D0*A4
 
672
      B5     = 5.D0*A5
 
673
      B6     = 6.D0*A6
 
674
      CNORM  = 10.D0/YSTP
 
675
 
 
676
C
 
677
C ... CHECK FOR OUTPUT START
 
678
C
 
679
      STARTL = YSTR - 0.5*YSTP
 
680
      FN     = 0.5D0
 
681
      XN     = IX_R8R8(FN,XSTR,XSTP)
 
682
      PN     = 10.D0*WPOLY(XN)
 
683
 
 
684
      DELT   = PN - STARTL
 
685
      IF (DELT.GT.0.) THEN
 
686
          M1     = DELT/YSTP + 1.
 
687
          TMP1   = PN + YSTP
 
688
          TM     = PN
 
689
      ELSE
 
690
          M1     = 1
 
691
          TMP1   = STARTL + YSTP
 
692
          TM     = STARTL
 
693
      END IF
 
694
C
 
695
C ... CHECK FOR INPUT START
 
696
C
 
697
      DO 10 N = 1,NX
 
698
          FN     = FN + 1.D0
 
699
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
 
700
          PNP1   = 10.D0*WPOLY(XNP1)
 
701
 
 
702
          IF (PNP1.GT.TM) GO TO 20
 
703
          PN     = PNP1
 
704
   10 CONTINUE
 
705
      N      = NX + 1
 
706
      GO TO 30
 
707
   20 C      = X(N)
 
708
   30 VAL    = 0.
 
709
 
 
710
      COVX   = XNP1 - XN
 
711
      COV    = (PNP1 - PN)/COVX
 
712
C
 
713
C ... REBINNING LOOP
 
714
C
 
715
      MPT    = 1
 
716
      MST    = 1
 
717
      K      = 0
 
718
   35 CONTINUE
 
719
          IF ( .NOT. (MPT.LE.NY)) GO TO 80
 
720
          K      = K + 1
 
721
          IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60
 
722
C 
 
723
   40     IF (PNP1.EQ.PN) GO TO 90
 
724
          VAL    = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV
 
725
          IF (TMP1.LT.PNP1) GO TO 50
 
726
          N      = N + 1
 
727
          IF (N.GT.NX) GO TO 60
 
728
 
 
729
          PN     = PNP1
 
730
          XN     = XNP1
 
731
 
 
732
          FN     = FN + 1.D0
 
733
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
 
734
          PNP1   = 10.D0*WPOLY(XNP1)
 
735
 
 
736
          C      = X(N)
 
737
          COV    = (PNP1 - PN)/(XNP1-XN)
 
738
 
 
739
          GO TO 40
 
740
   50     TM     = TMP1
 
741
          TMP1   = TMP1 + YSTP
 
742
   60     FACTOR = CNORM*DPOLY(IX_R8R8(FN,XSTR,XSTP))
 
743
C          FACTOR = FACTOR*10.D0*DABS(DPOLY(IX_R8R8(FN,XSTR,XSTP)))/YSTP
 
744
          Y(MPT) = VAL*FACTOR
 
745
          MPT    = MPT + 1
 
746
          VAL    = 0.
 
747
          GOTO 35
 
748
   80 CONTINUE
 
749
   90 CONTINUE
 
750
      DO 100 K = MPT + 1,NY
 
751
          Y(K)   = 0.
 
752
  100 CONTINUE
 
753
      RETURN
 
754
      END
 
755
 
 
756
      SUBROUTINE LREBIN(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP)
 
757
C
 
758
C
 
759
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
760
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
761
C                                         all rights reserved
 
762
C
 
763
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  19:08 - 6 DEC 1987
 
764
C
 
765
C.LANGUAGE: F77+ESOext
 
766
C
 
767
C.AUTHOR: D.PONZ
 
768
C
 
769
C.KEYWORDS:
 
770
C  SPECTROSCOPY, CASPEC, REBINNING
 
771
C
 
772
C.PURPOSE:
 
773
C  DATA REBINNING
 
774
C  CONVERSION FROM SAMPLES INTO WAVELENGTHS
 
775
C
 
776
C.ALGORITHM:
 
777
C  LINEAR INTERPOLATION. THE TRANSFORMATION WAVELENGTH INTO
 
778
C  SAMPLE NUMBER IS DEFINED BY
 
779
C  S = A0+A1*T+A2*T**2
 
780
C
 
781
C.INPUT/OUTPUT
 
782
C  INPUT ARGUMENTS
 
783
C
 
784
C  X(NX) REAL*4 INPUT ARRAY
 
785
C  NX  INTG*4 NO OF SAMPLES
 
786
C  XSTR  REAL*4 STARTING VALUE OF INPUT
 
787
C  XSTP  REAL*4 INPUT STEP
 
788
C  AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION
 
789
C  NC  REAL*4 NO. OF COEFFICIENTS
 
790
C  NY  INTG*4 NO. OF OUTPUT VALUES
 
791
C  YSTR  REAL*4 STARTING VALUE OF THE INPUT
 
792
C  YSTP  REAL*4 STEP OF THE OUTPUT
 
793
C
 
794
C  OUTPUT ARGUMENTS
 
795
C
 
796
C  Y(NY) REAL*4 OUTPUT ARRAY
 
797
C
 
798
C-------------------------------------------------------------
 
799
C
 
800
C
 
801
      IMPLICIT NONE
 
802
      INTEGER NX, NY, NC, N, M1
 
803
      REAL X(NX),Y(NY)
 
804
      DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,FN
 
805
      DOUBLE PRECISION FACTOR,CNORM, START, PN, DELT, STARTL
 
806
      DOUBLE PRECISION TM, TMP1
 
807
      DOUBLE PRECISION IX_R8R8
 
808
      DOUBLE PRECISION XSTR, XSTP, YSTR, YSTP
 
809
C
 
810
      INTEGER K, MPT, MST
 
811
      DOUBLE PRECISION PNP1, C, VAL, COV, DSQRT, DMIN1, DMAX1
 
812
      DOUBLE PRECISION XN, XNP1, COVX
 
813
C
 
814
C WPOLY(FL) = FL*(A1+FL*(A2+FL*(A3+FL*(A4+FL*A5))))
 
815
C
 
816
C ... ASSIGN COEFFS
 
817
C
 
818
      A0     = AC(1)
 
819
      A1     = -AC(2)
 
820
      A2     = AC(3)
 
821
      A3     = A1*A1
 
822
      A4     = 4.D0*A2
 
823
      A5     = 2.D0*A2
 
824
      CNORM  = 10.D0/ (A5*YSTP)
 
825
C A0 = AC(1)
 
826
C A1 = AC(2)
 
827
C A2 = AC(3)
 
828
C
 
829
C ... CHECK FOR OUTPUT START
 
830
C
 
831
      START  = YSTR - 0.5*YSTP
 
832
C STARTL= START - A0
 
833
      STARTL = START
 
834
      M1     = 1
 
835
      FN     = 0.5D0
 
836
      XN     = IX_R8R8(FN,XSTR,XSTP)
 
837
 
 
838
C PN    = WPOLY(IX_R8R8(FN,XSTR,XSTP))
 
839
      PN     = 10.D0* (A1+DSQRT(A3-A4* (A0-XN)))/A5
 
840
      DELT   = PN - STARTL
 
841
      IF (DELT.GT.0.) THEN
 
842
          M1     = DELT/YSTP + 1.
 
843
          TMP1   = STARTL + M1*YSTP
 
844
          TM     = TMP1 - YSTP
 
845
      ELSE
 
846
          TMP1   = STARTL + YSTP
 
847
          TM     = STARTL
 
848
      END IF
 
849
C
 
850
C ... CHECK FOR INPUT START
 
851
C
 
852
      DO 10 N = 1,NX
 
853
          FN     = FN + 1.D0
 
854
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
 
855
 
 
856
C   PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP))
 
857
          PNP1   = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5
 
858
          IF (PNP1.GT.TM) GO TO 20
 
859
          PN     = PNP1
 
860
          XN     = XNP1
 
861
   10 CONTINUE
 
862
      N      = NX + 1
 
863
      GO TO 30
 
864
   20 C      = X(N)
 
865
   30 VAL    = 0.
 
866
      COVX   = XNP1 - XN
 
867
      COV    = (PNP1 - PN)/COVX
 
868
C
 
869
C ... REBINNING LOOP
 
870
C
 
871
      MPT    = 1
 
872
      MST    = 1
 
873
      K      = 0
 
874
C   35 CONTINUE
 
875
          IF ( .NOT. (MPT.LE.NY)) GO TO 80
 
876
          K      = K + 1
 
877
          IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60
 
878
   40     VAL    = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV
 
879
          IF (TMP1.LT.PNP1) GO TO 50
 
880
          N      = N + 1
 
881
          IF (N.GT.NX) GO TO 60
 
882
          PN     = PNP1
 
883
          XN     = XNP1
 
884
          FN     = FN + 1.D0
 
885
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
 
886
C         PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP))
 
887
          PNP1   = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5
 
888
 
 
889
          C      = X(N)
 
890
 
 
891
          COVX   = XNP1 - XN
 
892
          COV    = (PNP1 - PN)/COVX
 
893
 
 
894
          GO TO 40
 
895
   50     TM     = TMP1
 
896
          TMP1   = TMP1 + YSTP
 
897
   60     FACTOR = CNORM*(DSQRT(A3-A4* (A0-XNP1))-
 
898
     +             DSQRT(A3-A4* (A0-XNP1-1.D0)))
 
899
          Y(MPT) = VAL*FACTOR
 
900
          MPT    = MPT + 1
 
901
          VAL    = 0.
 
902
   80 CONTINUE
 
903
      DO 90 K = MPT + 1,NY
 
904
          Y(K)   = 0.
 
905
   90 CONTINUE
 
906
      RETURN
 
907
      END
 
908