← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necmerge.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/necmerge.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2008 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C 
 
30
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
31
C                                         all rights reserved
 
32
C
 
33
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  22:26 - 3 DEC 1987
 
34
C
 
35
C.LANGUAGE: F77+ESOext
 
36
C
 
37
C.AUTHOR: D.PONZ
 
38
C
 
39
C
 
40
C.IDENTIFICATION
 
41
C
 
42
C  program  ECHMERG 
 
43
C
 
44
C.MODIFICATIONS
 
45
C
 
46
C  910128   P. Ballester   Define and set variables for mmake -i
 
47
C  991005   S.Wolf         ECHMR2O added: does an optimal merging.
 
48
C
 
49
C.KEYWORDS
 
50
C
 
51
C  ECHELLE, CASPEC, ORDER MERGING
 
52
C
 
53
C.PURPOSE
 
54
C
 
55
C  OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM
 
56
C  SAMPLED IN THE SPACE WAVELENGTH-ORDER.
 
57
C
 
58
C.ALGORITHM
 
59
C
 
60
C  OVERLAPPING REGIONS OF CONSECUTIVE ORDERS ARE PROCESSED IN TWO WAYS:
 
61
C  - CONCATENATION OF CONSECUTIVE ORDERS, USING THE MIDDLE POINT AS
 
62
C    CONCATENATING POSITION (METHOD='CONCATENATE', DEFAULT)
 
63
C  - AVERAGE OF OVERLAPPING REGION (METHOD='AVERAGE')
 
64
C  - WEIGHTED AVERAGE OF OVERLAPPING REGION (METHOD='OPTIMAL')
 
65
C    spmerged = (w1*sp1 + w2*sp2)/(w1+w2)
 
66
C  - NO MERGING. INDIVIDUAL ORDERS ARE WRITTEN IN DIFFERENT FILES
 
67
C    (METHOD = 'NOMERGE')
 
68
C  - WEIGHTED AVERAGE. WEIGHTS PROPORTIONAL TO SINC**2 (METHOD='SINC')
 
69
C  THE METHOD 'NOMERGE' WILL PRODUCE AS MANY FILES AS DEFINED ORDERS
 
70
C   WITH NAMES xxxyyyy, xxxx IS THE OUTPUT FILE NAME, yyyy IS THE 
 
71
C   ORDER NUMBER
 
72
C
 
73
C.INPUT/OUTPUT
 
74
C
 
75
C  MERGE/ECHELLE INPUT OUTPUT  [CONC]
 
76
C  MERGE/ECHELLE INPUT OUTPUT  ORD1,ORD2 NOCONCAT
 
77
C  MERGE/ECHELLE INPUT OUTPUT  DELTA     AVERAGE
 
78
C  MERGE/ECHELLE INPUT OUTPUT  DELTA     OPTIMAL WEIGHTS
 
79
C  or
 
80
C  MERGE/ECHELLE INPUT OUTPUT  TABLE SINC
 
81
C
 
82
C 081030        last modif
 
83
C 
 
84
C------------------------------------------------------------------------
 
85
C 
 
86
      PROGRAM ECHMRG
 
87
C 
 
88
      IMPLICIT  NONE
 
89
C
 
90
      INTEGER   MADRID
 
91
      INTEGER   I, II1, II2, IORD1
 
92
      INTEGER   IORD2
 
93
      INTEGER   ISTAT
 
94
      INTEGER   NAXISA, NAXISB, NAXISW
 
95
      INTEGER NPIXA(3),NPIXB(3),NPIXW(3),NPTOT(100),ICOL(3),KUN,KNUL
 
96
      INTEGER INIMA, OUTIMA, VARIMA, WGTIMA, TID
 
97
C 
 
98
      INTEGER*8 IPNTRA, IPNTRB, IPNTRV, IPNTRW              !pointers ...
 
99
C   
 
100
      DOUBLE PRECISION  DEL, WEND, WINIT, WSTEP
 
101
      DOUBLE PRECISION  WSTART(100)
 
102
      DOUBLE PRECISION  STARTA(3),STEPA(3)
 
103
      DOUBLE PRECISION  STARTB(3),STEPB(3)
 
104
      DOUBLE PRECISION  STARTW(3),STEPW(3)
 
105
C
 
106
      REAL CUT(4),XMIN,XMAX,WRANGE(2),VAL(3),K(100),A(100)
 
107
      REAL RO(100)
 
108
C
 
109
      CHARACTER*60 INFRM, OUTFRM, VARFRM, WGTFRM
 
110
      CHARACTER OUTFIL*12,TABLE*12
 
111
      CHARACTER METHOD*1,WS*5
 
112
      CHARACTER IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64
 
113
      CHARACTER IDENTW*72,CUNITW*64
 
114
      CHARACTER*16 LABEL(3)
 
115
C
 
116
      LOGICAL NULL(3)
 
117
C
 
118
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
119
      COMMON /VMR/MADRID(1)
 
120
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
121
C
 
122
      DATA CUNITB/'FLUX            WAVELENGTH'/
 
123
      DATA LABEL(1)/'ORDER'/,LABEL(2)/'KFIT'/,LABEL(3)/'AFIT'/
 
124
C
 
125
C ... INITIALIZE SYSTEM
 
126
C
 
127
      CALL STSPRO('ECHMRG')
 
128
      INFRM = ' '
 
129
      OUTFRM = ' '
 
130
      VARFRM = ' '
 
131
      WGTFRM = ' '
 
132
      CALL STKRDC('P1',1,1,60,I,INFRM,KUN,KNUL,ISTAT)
 
133
      CALL STKRDC('P2',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT)
 
134
      CALL STKRDC('P6',1,1,60,I,VARFRM,KUN,KNUL,ISTAT)
 
135
      CALL STKRDC('P5',1,1,60,I,WGTFRM,KUN,KNUL,ISTAT)
 
136
      CALL CLNFRA(INFRM,INFRM,0)
 
137
      CALL CLNFRA(OUTFRM,OUTFRM,0)
 
138
      CALL CLNFRA(VARFRM,VARFRM,0)
 
139
      CALL CLNFRA(WGTFRM,WGTFRM,0)
 
140
      CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT)
 
141
      CALL FORUPC(METHOD,METHOD)
 
142
      IF (METHOD.EQ.'S') 
 
143
     .    CALL STKRDC('P3',1,1,8,I,TABLE,KUN,KNUL,ISTAT)
 
144
      CALL STKRDR('INPUTR',1,2,I,WRANGE,KUN,KNUL,ISTAT)
 
145
      IF (METHOD.EQ.'A' .OR. METHOD.EQ.'O') THEN
 
146
          WINIT  = 0.
 
147
          WEND   = WINIT
 
148
          DEL    = WRANGE(1)
 
149
      ELSE
 
150
          WINIT  = WRANGE(1)
 
151
          WEND   = WRANGE(2)
 
152
      END IF
 
153
C
 
154
C ... MAP INPUT FRM
 
155
C
 
156
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
157
     .            3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
 
158
     .            IPNTRA,INIMA,ISTAT)
 
159
      CALL STDRDD(INIMA,'WSTART',1,NPIXA(2),I,
 
160
     .            WSTART,KUN,KNUL,ISTAT)
 
161
      CALL STDRDI(INIMA,'NPTOT',1,NPIXA(2),I,NPTOT,
 
162
     .            KUN,KNUL,ISTAT)
 
163
C ... CHECK IF NPTOT(NPIXA(2)) IS 0 
 
164
      IF (NPTOT(NPIXA(2)).EQ.0) THEN  
 
165
          NPIXA(2) = NPIXA(2) - 1
 
166
      END IF
 
167
      WSTEP  = STEPA(1)
 
168
C
 
169
C ... MAP WEIGHT FRM
 
170
C
 
171
      IF (METHOD.EQ.'O') THEN
 
172
         CALL STIGET(WGTFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
173
     .        3,NAXISW,NPIXW,STARTW,STEPW,IDENTW,CUNITW,
 
174
     .        IPNTRW,WGTIMA,ISTAT)
 
175
         IF (NAXISW.NE.NAXISA.OR.
 
176
     .        STARTW(1).NE.STARTA(1).OR.STARTW(2).NE.STARTA(2).OR.
 
177
     .        STEPW(1).NE.STEPA(1).OR.STEPW(1).NE.STEPA(1)) THEN
 
178
            GOTO 9998
 
179
         END IF
 
180
      END IF
 
181
      IF (METHOD.EQ.'N') THEN
 
182
C
 
183
C ... NO MERGE - PRODUCE ONE FILE PER ORDER
 
184
C
 
185
          IF (WINIT.EQ.WEND .AND. WINIT.LT.0.5) THEN
 
186
              IORD1  = 1
 
187
              IORD2  = NPIXA(2)
 
188
          ELSE
 
189
              IORD1  = MAX(NINT(WINIT),1)
 
190
              IORD2  = MIN(NINT(WEND),NPIXA(2))
 
191
          END IF
 
192
          OUTFRM = OUTFRM(1:59)//' '
 
193
          II1    = INDEX(OUTFRM,' ') - 1
 
194
          II1    = MIN(II1,4)
 
195
          OUTFIL = OUTFRM(1:II1)
 
196
          DO 10 I = IORD1,IORD2
 
197
              II2    = 10000 + I
 
198
              WRITE (WS,9000) II2
 
199
              OUTFIL(II1+1:II1+4) = WS(2:5)
 
200
              NAXISB    = 1
 
201
              NPIXB(1)  = NPTOT(I)
 
202
              NPIXB(2)  = 1
 
203
              STARTB(1) = WSTART(I)
 
204
              STARTB(2) = 1.
 
205
              STEPB(1)  = WSTEP
 
206
              STEPB(2)  = 0.
 
207
              WRITE (IDENTB,9010) I
 
208
              IDENTB(11:72) = IDENTA(1:62)
 
209
              CALL STIPUT(OUTFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
210
     .                    NAXISB,NPIXB,STARTB,STEPB,
 
211
     .                     IDENTB,CUNITB,IPNTRB,OUTIMA,ISTAT)
 
212
              CALL COPY(MADRID(IPNTRA),NPIXA(1),NPIXA(2),MADRID(IPNTRB),
 
213
     .                  NPIXB(1),I,XMIN,XMAX)
 
214
              CUT(1) = XMIN
 
215
              CUT(2) = XMAX
 
216
              CUT(3) = XMIN
 
217
              CUT(4) = XMAX
 
218
              CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
 
219
              CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)
 
220
              CALL STFCLO(OUTIMA,ISTAT)
 
221
              CALL STTPUT('File '//OUTFIL//' created ...',ISTAT)
 
222
   10     CONTINUE
 
223
 
 
224
      ELSE
 
225
C
 
226
C ... MAP OUTPUT FRAME
 
227
C
 
228
          IF (WINIT.EQ.WEND) THEN
 
229
              WINIT  = WSTART(1)
 
230
              WEND   = WSTART(NPIXA(2)) + (NPTOT(NPIXA(2))-1)*WSTEP
 
231
          END IF
 
232
          NAXISB    = 1
 
233
          NPIXB(1)  = (WEND-WINIT)/WSTEP + 1
 
234
          NPIXB(2)  = 1
 
235
          STARTB(1) = WINIT
 
236
          STARTB(2) = 1.
 
237
          STEPB(1)  = WSTEP
 
238
          STEPB(2)  = 0.
 
239
          CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
240
     .                    NAXISB,NPIXB,STARTB,STEPB,
 
241
     .                    IDENTA,CUNITB,IPNTRB,OUTIMA,ISTAT)
 
242
C
 
243
C ... METHODS OF OVERLAPPING
 
244
C
 
245
          IF (METHOD.EQ.'C') CALL ECHMR1(MADRID(IPNTRA),
 
246
     +                            NPIXA(1),NPIXA(2),
 
247
     +                            STEPA,WSTART,NPTOT,MADRID(IPNTRB),
 
248
     +                            NPIXB,STARTB,XMIN,XMAX)
 
249
          IF (METHOD.EQ.'A') CALL ECHMR2(MADRID(IPNTRA),
 
250
     +                            NPIXA(1),NPIXA(2),
 
251
     +                            STEPA,WSTART,NPTOT,MADRID(IPNTRB),
 
252
     +                            NPIXB,STARTB,XMIN,XMAX,DEL)
 
253
          IF (METHOD.EQ.'O') THEN
 
254
             CALL STIPUT(VARFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
255
     .            NAXISB,NPIXB,STARTB,STEPB,
 
256
     .            IDENTA,CUNITB,IPNTRV,VARIMA,ISTAT)
 
257
             CALL ECHMR2O(MADRID(IPNTRA),MADRID(IPNTRW),
 
258
     +            NPIXA(1),NPIXA(2),
 
259
     +            STEPA,WSTART,NPTOT,MADRID(IPNTRB),MADRID(IPNTRV),
 
260
     +            NPIXB,STARTB,XMIN,XMAX,DEL)
 
261
             CALL DSCUPT(INIMA,VARIMA,' ',ISTAT)
 
262
             
 
263
          ENDIF
 
264
          IF (METHOD.EQ.'S') THEN
 
265
C
 
266
C ... READ RIPPLE TABLE
 
267
C
 
268
              CALL TBTOPN(TABLE,F_IO_MODE,TID,ISTAT)
 
269
              CALL TBLSER(TID,LABEL(1),ICOL(1),ISTAT)
 
270
              CALL TBLSER(TID,LABEL(2),ICOL(2),ISTAT)
 
271
              CALL TBLSER(TID,LABEL(3),ICOL(3),ISTAT)
 
272
              DO 20 I = 1,NPIXA(2)
 
273
                  CALL TBRRDR(TID,I,3,ICOL,VAL,NULL,ISTAT)
 
274
                  RO(I)  = VAL(1)
 
275
                  K(I)   = VAL(2)
 
276
                  A(I)   = VAL(3)
 
277
   20         CONTINUE
 
278
              CALL ECHMR3(MADRID(IPNTRA),
 
279
     +                      NPIXA(1),NPIXA(2),STEPA,WSTART,
 
280
     +                      NPTOT,MADRID(IPNTRB),NPIXB,STARTB,XMIN,XMAX,
 
281
     +                      RO,K,A)
 
282
          END IF
 
283
C
 
284
C ... WRITE PROCESS DESCRIPTORS
 
285
C
 
286
          CUT(1) = XMIN
 
287
          CUT(2) = XMAX
 
288
          CUT(3) = XMIN
 
289
          CUT(4) = XMAX
 
290
          CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
 
291
          CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)
 
292
      END IF
 
293
      GOTO 9999
 
294
 9998 CALL STTPUT('Error: Bad weight map! ',ISTAT)
 
295
      CALL STTPUT
 
296
     + ('(START, STEP and NAXIS same as in INPUT frame !?)',ISTAT)
 
297
 9999 CALL STSEPI
 
298
 9000 FORMAT (I5)
 
299
 9010 FORMAT ('ORDER:',I3)
 
300
      END
 
301
 
 
302
      SUBROUTINE ECHMR1(X,NPIXA1,NPIXA2,STEPA,WI,NP,Y,NY,YSTR,
 
303
     +                  XMIN,XMAX)
 
304
C
 
305
C MERGE THE ORDERS, USING SIMPLE CONCATENATION IN THE MIDDLE POINT
 
306
      IMPLICIT NONE
 
307
C
 
308
      INTEGER  NY,J2,I,J1,JOFF,J,JJ
 
309
      INTEGER  NPIXA1,NPIXA2
 
310
      INTEGER  NP(NPIXA2)
 
311
C 
 
312
      DOUBLE PRECISION STEPA(2),WI(NPIXA2),YSTR
 
313
C 
 
314
      REAL     X(NPIXA1,NPIXA2)
 
315
      REAL     XMIN,XMAX,RVAL,WEN1,WST2,WSTART
 
316
      REAL     Y(NY)
 
317
C 
 
318
      DOUBLE PRECISION WEND, YSTEP, YEND, WSTR
 
319
C
 
320
      XMIN   = 0.
 
321
      XMAX   = 0.
 
322
      DO 10 I = 1,NY
 
323
          Y(I)   = 0.
 
324
   10 CONTINUE
 
325
      YSTEP  = STEPA(1)
 
326
      YEND   = YSTR + (NY-1)*YSTEP
 
327
C
 
328
C ... ITERATION ON ORDERS
 
329
C
 
330
      WEND   = 0.
 
331
      DO 30 I = 1,NPIXA2
 
332
          WSTR   = DMAX1(WI(I),WEND+YSTEP)
 
333
          IF (I.EQ.NPIXA2) THEN
 
334
              WEND   = WI(I) + (NP(I)-1)*YSTEP
 
335
          ELSE
 
336
              WEN1   = WI(I) + (NP(I)-1)*YSTEP
 
337
              WST2   = WI(I+1)
 
338
              IF (WST2.LT.WEN1) THEN
 
339
                  WEND   = 0.5* (WEN1+WST2)
 
340
              ELSE
 
341
                  WEND   = WEN1
 
342
              END IF
 
343
          END IF
 
344
 
 
345
          IF (WSTR.GE.YEND) RETURN
 
346
          IF (WEND.GT.YSTR) THEN
 
347
C
 
348
C ... ITERATION ON WAVELENGTHS
 
349
C
 
350
              WSTART = DMAX1(YSTR,WSTR)
 
351
              WEND   = DMIN1(WEND,YEND)
 
352
              J1     = NINT((WSTART-WI(I))/YSTEP) + 1
 
353
              J2     = NINT((WEND-WI(I))/YSTEP) + 1
 
354
              JOFF   = NINT((WI(I)-YSTR)/YSTEP)
 
355
              DO 20 J = J1,J2
 
356
                  JJ     = J + JOFF
 
357
                  IF (JJ.GT.0) THEN
 
358
                      RVAL = X(J,I)
 
359
                      IF (RVAL .GT. XMAX) XMAX = RVAL
 
360
                      IF (RVAL .LT. XMIN) XMIN = RVAL
 
361
                      Y(JJ) = RVAL
 
362
                  END IF
 
363
 
 
364
   20         CONTINUE
 
365
          END IF
 
366
 
 
367
   30 CONTINUE
 
368
      RETURN
 
369
 
 
370
      END
 
371
 
 
372
      SUBROUTINE ECHMR2(X,NPIXA1,NPIXA2,STEPA,WI,NP,Y,NY,YSTR,
 
373
     +                  XMIN,XMAX,DEL)
 
374
C
 
375
C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION
 
376
C
 
377
      IMPLICIT NONE
 
378
C
 
379
      INTEGER NPIXA1, NPIXA2, NY
 
380
      INTEGER NP(NPIXA2)
 
381
      INTEGER IORD1, IORD2, IPIX, IPIX1, IPIX2
 
382
C 
 
383
      REAL X(NPIXA1,NPIXA2)
 
384
      REAL Y(NY),XMIN,XMAX,RVAL
 
385
C 
 
386
      DOUBLE PRECISION STEPA(2),YSTR,DEL,WI(NPIXA2)
 
387
      DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2
 
388
C
 
389
C
 
390
      XMIN   = 0.
 
391
      XMAX   = 0.
 
392
      YSTEP  = STEPA(1)
 
393
      IORD1  = 1
 
394
      IORD2  = 2
 
395
      W0     = WI(IORD2) + DEL
 
396
      W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
 
397
C
 
398
C ... ITERATION ON ORDERS
 
399
C
 
400
      DO 20 IPIX = 1,NY
 
401
          Y(IPIX) = 0.0
 
402
          WL     = YSTR + (IPIX-1)*YSTEP
 
403
          IF (WL.LT.W0) THEN
 
404
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
405
              RVAL = X(IPIX1,IORD1)
 
406
              IF (RVAL .GT. XMAX) XMAX = RVAL
 
407
              IF (RVAL .LT. XMIN) XMIN = RVAL
 
408
              Y(IPIX) = RVAL
 
409
          ELSE IF (WL.LT.W1) THEN
 
410
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
411
              IPIX2  = NINT((WL-WI(IORD2))/YSTEP) + 1
 
412
              P2     = (WL-W0)/ (W1-W0)
 
413
              P1     = 1.D0 - P2
 
414
              IF (X(IPIX1,IORD1).LE.0.0) THEN
 
415
                  P2     = 1.D0
 
416
                  P1     = 0.D0
 
417
              END IF
 
418
              IF (X(IPIX2,IORD2).LE.0.0) THEN
 
419
                  P2     = 0.D0
 
420
                  P1     = 1.D0
 
421
              END IF
 
422
              RVAL = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2
 
423
              IF (RVAL .GT. XMAX) XMAX = RVAL
 
424
              IF (RVAL .LT. XMIN) XMIN = RVAL
 
425
              Y(IPIX) = RVAL
 
426
          ELSE
 
427
              IORD1  = IORD1 + 1
 
428
              IF (IORD1.GT.NPIXA2) RETURN
 
429
              IORD2  = IORD2 + 1
 
430
              IF (IORD2.GT.NPIXA2) THEN
 
431
                  W0     = 1.E20
 
432
              ELSE
 
433
                  W0     = WI(IORD2) + DEL
 
434
              END IF
 
435
              W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
 
436
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
437
              RVAL = X(IPIX1,IORD1)
 
438
              IF (RVAL .GT. XMAX) XMAX = RVAL
 
439
              IF (RVAL .LT. XMIN) XMIN = RVAL
 
440
              Y(IPIX) = RVAL
 
441
          END IF
 
442
   20 CONTINUE
 
443
      RETURN
 
444
      END
 
445
 
 
446
      SUBROUTINE ECHMR2O(X,WGT,NPIXA1,NPIXA2,STEPA,WI,NP,Y,VAR,NY
 
447
     $     ,YSTR,XMIN,XMAX,DEL)
 
448
C
 
449
C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION USING WEIGHTS
 
450
C
 
451
      IMPLICIT NONE
 
452
C
 
453
      INTEGER NPIXA1, NPIXA2, NY
 
454
      INTEGER NP(NPIXA2)
 
455
      REAL X(NPIXA1,NPIXA2), WGT(NPIXA1,NPIXA2)
 
456
      DOUBLE PRECISION STEPA(2),YSTR,DEL,WI(NPIXA2)
 
457
      REAL Y(NY), VAR(NY),XMIN,XMAX,RVAL
 
458
C
 
459
      INTEGER IORD1, IORD2, IPIX, IPIX1, IPIX2, IBAD
 
460
      DOUBLE PRECISION YSTEP, WG1, WG2, W0, W1, WL
 
461
      CHARACTER*80 MES
 
462
C
 
463
      IBAD   = 0
 
464
      XMIN   = 0.
 
465
      XMAX   = 0.
 
466
      YSTEP  = STEPA(1)
 
467
      IORD1  = 1
 
468
      IORD2  = 2
 
469
      W0     = WI(IORD2) + DEL
 
470
      W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
 
471
C
 
472
C ... ITERATION ON ORDERS
 
473
C
 
474
      DO 20 IPIX = 1,NY
 
475
          Y(IPIX) = 0.0
 
476
          VAR(IPIX) = 0.0
 
477
          WL     = YSTR + (IPIX-1)*YSTEP
 
478
          IF (WL.LT.W0) THEN
 
479
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
480
              RVAL = X(IPIX1,IORD1)
 
481
              IF (RVAL .GT. XMAX) XMAX = RVAL
 
482
              IF (RVAL .LT. XMIN) XMIN = RVAL
 
483
              Y(IPIX) = RVAL
 
484
              VAR(IPIX) = WGT(IPIX1,IORD1)
 
485
              IF (VAR(IPIX).NE.0.0) VAR(IPIX) = 1.0/VAR(IPIX)
 
486
          ELSE IF (WL.LT.W1) THEN
 
487
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
488
              IPIX2  = NINT((WL-WI(IORD2))/YSTEP) + 1
 
489
              WG1 = WGT(IPIX1,IORD1)
 
490
              WG2 = WGT(IPIX2,IORD2)
 
491
 
 
492
              IF (WG1.LT.1e-10 .AND. WG2.LT.1e-10) THEN
 
493
                 Y(IPIX) = 0.0
 
494
                 VAR(IPIX) = 0.0
 
495
                 IBAD = IBAD + 1
 
496
              ELSE
 
497
                 VAR(IPIX) = 1.0/(WG1 + WG2)
 
498
                 Y(IPIX) = X(IPIX1,IORD1)*WG1 + X(IPIX2,IORD2)*WG2
 
499
                 RVAL = Y(IPIX) * VAR(IPIX)
 
500
                 IF (RVAL .GT. XMAX) XMAX = RVAL
 
501
                 IF (RVAL .LT. XMIN) XMIN = RVAL
 
502
                 Y(IPIX) = RVAL
 
503
              ENDIF
 
504
          ELSE
 
505
              IORD1  = IORD1 + 1
 
506
              IF (IORD1.GT.NPIXA2) GOTO 99
 
507
              IORD2  = IORD2 + 1
 
508
              IF (IORD2.GT.NPIXA2) THEN
 
509
                  W0 = 1.E20
 
510
              ELSE
 
511
                  W0 = WI(IORD2) + DEL
 
512
              END IF
 
513
              W1 = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
 
514
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
 
515
              RVAL = X(IPIX1,IORD1)
 
516
              IF (RVAL .GT. XMAX) XMAX = RVAL
 
517
              IF (RVAL .LT. XMIN) XMIN = RVAL
 
518
              Y(IPIX) = RVAL
 
519
              VAR(IPIX) = WGT(IPIX1,IORD1)
 
520
              IF (VAR(IPIX).NE.0.0) VAR(IPIX) = 1.0/VAR(IPIX)
 
521
          END IF
 
522
 20   CONTINUE
 
523
 99   IF (IBAD.GT.0) THEN
 
524
         WRITE(MES,9099) IBAD
 
525
         CALL STTPUT(MES,IBAD)
 
526
      END IF      
 
527
      RETURN
 
528
 
 
529
9099  FORMAT(I5,' undefined pixels ... set to 0.0!')
 
530
         
 
531
      END
 
532
 
 
533
      SUBROUTINE COPY(A,NA1,NA2,B,NB,I,XMIN,XMAX)
 
534
C
 
535
C COPY THE ORDER NUMBER I
 
536
C
 
537
      IMPLICIT NONE
 
538
C 
 
539
      INTEGER NA1, NA2, NB, I, J
 
540
C 
 
541
      REAL A(NA1,NA2),B(NB)
 
542
      REAL XMIN,XMAX,RVAL
 
543
C
 
544
      XMIN = 0.
 
545
      XMAX = 0.
 
546
C
 
547
      DO 10 J = 1,NB
 
548
          B(J)   = A(J,I)
 
549
          RVAL = B(J)
 
550
          IF (RVAL .GT. XMAX) XMAX = RVAL
 
551
          IF (RVAL .LT. XMIN) XMIN = RVAL
 
552
   10 CONTINUE
 
553
      RETURN
 
554
      END
 
555
 
 
556
      SUBROUTINE ECHMR3(X,NPIXA1,NPIXA2,STEPA,WI,NP,Y,NY,YSTR,
 
557
     +                    XMIN,XMAX,RORDER,K,A)
 
558
C
 
559
C      IMPLICIT NONE
 
560
C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION WITH
 
561
C WEIGHTS PROPORTIONAL TO THE SINC**2
 
562
C
 
563
      IMPLICIT NONE
 
564
C 
 
565
      INTEGER NY,IORD1,IPIX,NF,IORD,IPIX1,I,IFL,IORD2
 
566
      INTEGER NPIXA1,NPIXA2
 
567
      INTEGER NP(NPIXA2)
 
568
C 
 
569
      REAL    XMIN,XMAX,YSTEP,PI,DW,WL
 
570
      REAL    X(NPIXA1,NPIXA2)
 
571
      REAL    Y(NY),RORDER(NPIXA2)
 
572
      REAL    K(NPIXA2),A(NPIXA2)
 
573
C 
 
574
      DOUBLE PRECISION DK,DA,PA,DM,DC,DX,WEIGHT(3),FLUX(3),SW,SF
 
575
      DOUBLE PRECISION YSTR,STEPA(2),WI(NPIXA2)
 
576
C
 
577
      XMIN   = 0.
 
578
      XMAX   = 0.
 
579
      PI     = 0.     !  To be updated (ULTRIX installation)
 
580
      DW     = 0.     !  To be updated (ULTRIX installation)
 
581
      YSTEP  = STEPA(1)
 
582
      IORD1  = 1
 
583
      IORD2  = MIN(IORD1+2,NPIXA2)
 
584
C
 
585
C ... ITERATION ON OUTPUT WAVELENGTHS
 
586
C
 
587
      DO 30 IPIX = 1,NY
 
588
          Y(IPIX) = 0.
 
589
          FLUX(1) = 0.D0
 
590
          FLUX(2) = 0.D0
 
591
          FLUX(3) = 0.D0
 
592
          WEIGHT(1) = 0.D0
 
593
          WEIGHT(2) = 0.D0
 
594
          WEIGHT(3) = 0.D0
 
595
          NF     = 0
 
596
          WL     = YSTR + (IPIX-1)*YSTEP
 
597
C
 
598
C ... ITERATION ON ORDERS
 
599
C
 
600
          DO 10 IORD = IORD1,IORD2
 
601
              IPIX1  = (WL-WI(IORD))/YSTEP + 1
 
602
              IF (IPIX1.GT.5 .AND. IPIX1.LT. (NP(IORD)-5)) THEN
 
603
C
 
604
C ... PIXEL IN THE ORDER RANGE
 
605
C
 
606
                  NF     = NF + 1
 
607
                  FLUX(NF) = X(IPIX1,IORD)
 
608
                  DK     = K(IORD)
 
609
                  DA     = A(IORD)
 
610
                  PA     = PI*DA
 
611
                  DM     = RORDER(IORD)
 
612
                  DC     = DM/DK
 
613
                  DX     = (PA*DM*DC)* (DW-1.D0/DC)
 
614
                  IF (DABS(DX).LT.1.D-10) THEN
 
615
                      WEIGHT(NF) = 1.D0
 
616
 
 
617
                  ELSE
 
618
                      WEIGHT(NF) = 1.D0/ (DSIN(DX)/DX)**4
 
619
                  END IF
 
620
 
 
621
              END IF
 
622
 
 
623
   10     CONTINUE
 
624
C
 
625
C ... AVERAGE OVER FLUXES
 
626
C
 
627
          IF (NF.GE.1) THEN
 
628
              SW     = 0.D0
 
629
              SF     = 0.D0
 
630
              DO 20 IFL = 1,NF
 
631
                  SW     = SW + WEIGHT(IFL)
 
632
                  SF     = SF + WEIGHT(IFL)*FLUX(IFL)
 
633
   20         CONTINUE
 
634
              Y(IPIX) = SF/SW
 
635
              XMIN   = MIN(XMIN,Y(IPIX))
 
636
              XMAX   = MAX(XMAX,Y(IPIX))
 
637
          END IF
 
638
C
 
639
C ... CHECK ON ORDER LIMITS
 
640
C
 
641
          IPIX1  = (WL-WI(IORD1))/YSTEP + 1
 
642
          IF (IPIX1.GE.NP(IORD1)-5) THEN
 
643
              IORD1  = IORD1 + 1
 
644
              IORD2  = MIN(IORD1+2,NPIXA2)
 
645
          END IF
 
646
 
 
647
          IF (IORD1.GT.NPIXA2) RETURN
 
648
   30 CONTINUE
 
649
      DO 40 I = IPIX + 1,NY
 
650
          Y(I)   = 0.
 
651
   40 CONTINUE
 
652
      RETURN
 
653
 
 
654
      END