C @(#)echmerg1.for 17.1.1.1 (ESO-IPG) 01/25/02 17:52:20 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM ECHMRG C + C .COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C .VERSION: 1.0 ESO-FORTRAN Conversion, AA 22:26 - 3 DEC 1987 C C .LANGUAGE: F77+ESOext C C .AUTHOR: D.PONZ C C C .IDENTIFICATION C C program ECHMERG C C .MODIFICATIONS C C 910128 P. Ballester Define and set variables for mmake -i C 981104 O. Stahl Modified for FEROS C C .KEYWORDS C C ECHELLE, CASPEC, ORDER MERGING C C .PURPOSE C C OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM C SAMPLED IN THE SPACE WAVELENGTH-ORDER. C C .ALGORITHM C C OVERLAPPING REGIONS OF CONSECUTIVE ORDERS ARE PROCESSED IN TWO WAYS: C - CONCATENATION OF CONSECUTIVE ORDERS, USING THE MIDDLE POINT AS C CONCATENATING POSITION (METHOD='CONCATENATE', DEFAULT) C - AVERAGE OF OVERLAPPING REGION (METHOD='AVERAGE') C - NO MERGING. INDIVIDUAL ORDERS ARE WRITTEN IN DIFFERENT FILES C (METHOD = 'NOMERGE') C - WEIGHTED AVERAGE. WEIGHTS PROPORTIONAL TO SINC**2 (METHOD='SINC') C THE METHOD 'NOMERGE' WILL PRODUCE AS MANY FILES AS DEFINED ORDERS C WITH NAMES xxxyyyy, xxxx IS THE OUTPUT FILE NAME, yyyy IS THE C ORDER NUMBER C C .INPUT/OUTPUT C C MERGE/ECHELLE INPUT OUTPUT [CONC] C MERGE/ECHELLE INPUT OUTPUT ORD1,ORD2 NOCONCAT C MERGE/ECHELLE INPUT OUTPUT DEL AVERAGE C OR C MERGE/ECHELLE INPUT OUTPUT TABLE SINC C C 010201 last modif C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C IMPLICIT NONE C INTEGER MADRID INTEGER I, II1, II2, IORD1 INTEGER IORD2 INTEGER*8 IPNTRA, IPNTRB INTEGER ISTAT INTEGER NAXISA, NAXISB C DOUBLE PRECISION DEL, WEND, WINIT, WSTEP C INTEGER NPIXA(3),NPIXB(3),NPTOT(100),ICOL(3),KUN,KNUL INTEGER INIMA, OUTIMA, TID DOUBLE PRECISION WSTART(100),STEPA(3) DOUBLE PRECISION STARTA(3),STEPB(3),STARTB(3) REAL CUT(4),XMIN,XMAX,WRANGE(2),VAL(3) REAL WAVE(100),ORDER(100),WAVE1(100),WAVE2(100) REAL WAVE3(100),WAVE4(100) CHARACTER*60 INFRM,OUTFRM CHARACTER OUTFIL*60 CHARACTER METHOD*1,WS*5 CHARACTER IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64 CHARACTER*16 LABEL(3) LOGICAL NULL(3) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C DATA CUNITB/'WAVELENGTH FLUX '/ DATA CUNITB/'FLUX WAVELENGTH'/ DATA LABEL(1)/'WAVE'/,LABEL(2)/'ORDER'/,LABEL(3)/'LNWAVE'/ C C ... INITIALIZE SYSTEM C CALL STSPRO('ECHMRG') INFRM = ' ' OUTFRM = ' ' CALL STKRDC('P1',1,1,60,I,INFRM,KUN,KNUL,ISTAT) CALL STKRDC('P2',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT) CALL CLNFRA(INFRM,INFRM,0) CALL CLNFRA(OUTFRM,OUTFRM,0) CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT) CALL FORUPC(METHOD,METHOD) CALL STKRDR('INPUTR',1,2,I,WRANGE,KUN,KNUL,ISTAT) C IF (METHOD.EQ.'A') THEN WINIT = 0. WEND = WINIT DEL = WRANGE(1) ELSE WINIT = WRANGE(1) WEND = WRANGE(2) END IF C IF (METHOD.EQ.'S') THEN WINIT = 0. WEND = WINIT DEL = WRANGE(1) END IF C C ... MAP INPUT FRM C CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . 3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, . IPNTRA,INIMA,ISTAT) CALL STDRDD(INIMA,'WSTART',1,NPIXA(2),I, . WSTART,KUN,KNUL,ISTAT) CALL STDRDI(INIMA,'NPTOT',1,NPIXA(2),I,NPTOT, . KUN,KNUL,ISTAT) C ... CHECK IF NPTOT(NPIXA(2)) IS 0 IF (NPTOT(NPIXA(2)).EQ.0) THEN NPIXA(2) = NPIXA(2) - 1 END IF WSTEP = STEPA(1) IF (METHOD.EQ.'N') THEN C C ... NO MERGE - PRODUCE ONE FILE PER ORDER C IF (WINIT.EQ.WEND .AND. WINIT.LT.0.5) THEN IORD1 = 1 IORD2 = NPIXA(2) ELSE IORD1 = MAX(NINT(WINIT),1) IORD2 = MIN(NINT(WEND),NPIXA(2)) END IF OUTFRM = OUTFRM(1:59)//' ' II1 = INDEX(OUTFRM,' ') - 1 II1 = MIN(II1,20) OUTFIL = OUTFRM(1:II1) DO 10 I = IORD1,IORD2 II2 = 10000 + I WRITE (WS,9000) II2 OUTFIL(II1+1:II1+4) = WS(2:5) NAXISB = 1 NPIXB(1) = NPTOT(I) NPIXB(2) = 1 STARTB(1) = WSTART(I) STARTB(2) = 1. STEPB(1) = WSTEP STEPB(2) = 0. WRITE (IDENTB,9010) I IDENTB(11:72) = IDENTA(1:62) CALL STIPUT(OUTFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NAXISB,NPIXB,STARTB,STEPB, . IDENTB,CUNITB,IPNTRB,OUTIMA,ISTAT) CALL COPY(MADRID(IPNTRA),NPIXA(1),NPIXA(2),MADRID(IPNTRB), . NPIXB(1),I,XMIN,XMAX) CUT(1) = XMIN CUT(2) = XMAX CUT(3) = XMIN CUT(4) = XMAX CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT) CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT) C CALL STFCLO(OUTIMA,ISTAT) CALL STTPUT('File '//OUTFIL//' created ...',ISTAT) 10 CONTINUE ELSE C C ... MAP OUTPUT FRAME C IF (WINIT.EQ.WEND) THEN WINIT = WSTART(1) WEND = WSTART(NPIXA(2)) + (NPTOT(NPIXA(2))-1)*WSTEP END IF NAXISB = 1 NPIXB(1) = (WEND-WINIT)/WSTEP + 1 NPIXB(2) = 1 STARTB(1) = WINIT STARTB(2) = 1. STEPB(1) = WSTEP STEPB(2) = 0. CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NAXISB,NPIXB,STARTB,STEPB, . IDENTA,CUNITB,IPNTRB,OUTIMA,ISTAT) C C ... METHODS OF OVERLAPPING C IF (METHOD.EQ.'C') CALL ECHMR1(MADRID(IPNTRA), + NPIXA(1),NPIXA(2),STARTA, + STEPA,WSTART,NPTOT,MADRID(IPNTRB), + NPIXB,STARTB,XMIN,XMAX) IF (METHOD.EQ.'A') CALL ECHMR2(MADRID(IPNTRA), + NPIXA(1),NPIXA(2),STARTA, + STEPA,WSTART,NPTOT,MADRID(IPNTRB), + NPIXB,STARTB,XMIN,XMAX,DEL) IF (METHOD.EQ.'S') THEN C C ... READ RIPPLE TABLE ------------------------------------------ C CALL TBTOPN('BLAZE',F_I_MODE,TID,ISTAT) CALL TBLSER(TID,LABEL(1),ICOL(1),ISTAT) CALL TBLSER(TID,LABEL(2),ICOL(2),ISTAT) CALL TBLSER(TID,LABEL(3),ICOL(3),ISTAT) C DO 20 I = 1,NPIXA(2)-1 CALL TBRRDR(TID,I,2,ICOL,VAL,NULL,ISTAT) IF(WSTART(1).GT.1000) THEN C looks like linear wavelength scale WAVE(I) = VAL(1) ELSE C looks like logarithmic wavelength scale WAVE(I) = VAL(3) ENDIF ORDER(I) = VAL(2) 20 CONTINUE DO 30 I=1,NPIXA(2) WAVE1(I) = WSTART(I)+DEL WAVE2(I) = WSTART(I)+(NPTOT(I)-1)*STEPA(1) - DEL 30 CONTINUE DO 60 I=1,NPIXA(2) WAVE3(I) = WAVE1(I) WAVE4(I) = WAVE2(I) 60 CONTINUE DO 40 I = 1,NPIXA(2)-2 WAVE3(I+1) = (WAVE(I)+WAVE(I+1))/2 - + (WAVE(I+1)-WAVE(I))/2*1.5 WAVE4(I+1) = (WAVE(I)+WAVE(I+1))/2 + + (WAVE(I+1)-WAVE(I))/2*1.5 40 CONTINUE DO 50 I=1,NPIXA(2) WAVE1(I) = AMAX1(WAVE1(I),WAVE3(I)) WAVE2(I) = AMIN1(WAVE2(I),WAVE4(I)) 50 CONTINUE CALL ECHMR3(MADRID(IPNTRA), + NPIXA(1),NPIXA(2),STARTA, + STEPA,WSTART,NPTOT,MADRID(IPNTRB), + NPIXB,STARTB,XMIN,XMAX,DEL,WAVE1,WAVE2) END IF C ... DONE WITH RIPPLE TABLE -------------------------------------- C C ... WRITE PROCESS DESCRIPTORS C CUT(1) = XMIN CUT(2) = XMAX CUT(3) = XMIN CUT(4) = XMAX CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT) CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT) END IF CALL STSEPI 9000 FORMAT (I5) 9010 FORMAT ('ORDER:',I3) END SUBROUTINE ECHMR1(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR, + XMIN,XMAX) C C MERGE THE ORDERS, USING SIMPLE CONCATENATION IN THE MIDDLE POINT IMPLICIT NONE C INTEGER NY,J2,I,J1,JOFF,J,JJ REAL XMIN,XMAX,WEN1,WST2,WSTART INTEGER NPIXA1,NPIXA2 INTEGER NP(NPIXA2) REAL X(NPIXA1,NPIXA2) DOUBLE PRECISION STARTA(2),STEPA(2),WI(NPIXA2),YSTR REAL Y(NY) DOUBLE PRECISION WEND, YSTEP, YEND, WSTR C XMIN = 0. XMAX = 0. DO 10 I = 1,NY Y(I) = 0. 10 CONTINUE YSTEP = STEPA(1) YEND = YSTR + (NY-1)*YSTEP C C ... ITERATION ON ORDERS C WEND = 0. DO 30 I = 1,NPIXA2 WSTR = DMAX1(WI(I),WEND+YSTEP) IF (I.EQ.NPIXA2) THEN WEND = WI(I) + (NP(I)-1)*YSTEP ELSE WEN1 = WI(I) + (NP(I)-1)*YSTEP WST2 = WI(I+1) IF (WST2.LT.WEN1) THEN WEND = 0.5* (WEN1+WST2) ELSE WEND = WEN1 END IF END IF IF (WSTR.GE.YEND) RETURN IF (WEND.GT.YSTR) THEN C C ... ITERATION ON WAVELENGTHS C WSTART = DMAX1(YSTR,WSTR) WEND = DMIN1(WEND,YEND) J1 = NINT((WSTART-WI(I))/YSTEP) + 1 J2 = NINT((WEND-WI(I))/YSTEP) + 1 JOFF = NINT((WI(I)-YSTR)/YSTEP) DO 20 J = J1,J2 JJ = J + JOFF IF (JJ.GT.0) THEN Y(JJ) = X(J,I) XMAX = AMAX1(XMAX,Y(JJ)) XMIN = AMIN1(XMIN,Y(JJ)) END IF 20 CONTINUE END IF 30 CONTINUE RETURN END SUBROUTINE ECHMR2(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR, + XMIN,XMAX,DEL) C C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION C IMPLICIT NONE C INTEGER NPIXA1, NPIXA2, NY INTEGER NP(NPIXA2) REAL X(NPIXA1,NPIXA2) DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2) REAL Y(NY), XMIN, XMAX C INTEGER I, IORD1, IORD2, IPIX, IPIX1, IPIX2 DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2, AMAX1, AMIN1 C XMIN = 0. XMAX = 0. DO 10 I = 1,NY Y(I) = 0. 10 CONTINUE YSTEP = STEPA(1) IORD1 = 1 IORD2 = 2 W0 = WI(IORD2) + DEL W1 = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL C C ... ITERATION ON ORDERS C DO 20 IPIX = 1,NY WL = YSTR + (IPIX-1)*YSTEP IF (WL.LT.W0) THEN IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 Y(IPIX) = X(IPIX1,IORD1) XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) ELSE IF (WL.LT.W1) THEN IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 IPIX2 = NINT((WL-WI(IORD2))/YSTEP) + 1 P2 = (WL-W0)/ (W1-W0) P1 = 1.D0 - P2 IF (X(IPIX1,IORD1).LE.0.0) THEN P2 = 1.D0 P1 = 0.D0 END IF IF (X(IPIX2,IORD2).LE.0.0) THEN P2 = 0.D0 P1 = 1.D0 END IF Y(IPIX) = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2 XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) ELSE IORD1 = IORD1 + 1 IF (IORD1.GT.NPIXA2) RETURN IORD2 = IORD2 + 1 IF (IORD2.GT.NPIXA2) THEN W0 = 1.E20 ELSE W0 = WI(IORD2) + DEL END IF W1 = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 Y(IPIX) = X(IPIX1,IORD1) XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) END IF 20 CONTINUE RETURN END SUBROUTINE COPY(A,NA1,NA2,B,NB,I,XMIN,XMAX) C C COPY THE ORDER NUMBER I C IMPLICIT NONE INTEGER NA1, NA2, NB, I, J REAL A(NA1,NA2),B(NB) REAL XMIN, XMAX C XMIN = 0. XMAX = 0. C DO 10 J = 1,NB B(J) = A(J,I) XMIN = AMIN1(B(J),XMIN) XMAX = AMAX1(B(J),XMAX) 10 CONTINUE RETURN END SUBROUTINE ECHMR3(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR, + XMIN,XMAX,DEL,WAVE1,WAVE2) C C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION C IMPLICIT NONE C INTEGER NPIXA1, NPIXA2, NY INTEGER NP(NPIXA2) REAL X(NPIXA1,NPIXA2) DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2) REAL Y(NY), XMIN, XMAX, WAVE1(NPIXA2), WAVE2(NPIXA2) C INTEGER I, IORD1, IORD2, IPIX, IPIX1, IPIX2 DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2, AMAX1, AMIN1 C XMIN = 0. XMAX = 0. DO 10 I = 1,NY Y(I) = 0. 10 CONTINUE YSTEP = STEPA(1) IORD1 = 1 IORD2 = 2 C W0 = WI(IORD2) + DEL C W1 = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL W0 = WAVE1(IORD2) W1 = WAVE2(IORD1) C C ... ITERATION ON ORDERS C DO 20 IPIX = 1,NY WL = YSTR + (IPIX-1)*YSTEP IF (WL.LT.W0) THEN IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 Y(IPIX) = X(IPIX1,IORD1) XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) ELSE IF (WL.LT.W1) THEN IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 IPIX2 = NINT((WL-WI(IORD2))/YSTEP) + 1 P2 = (WL-W0)/ (W1-W0) P1 = 1.D0 - P2 IF (X(IPIX1,IORD1).LE.0.0) THEN P2 = 1.D0 P1 = 0.D0 END IF IF (X(IPIX2,IORD2).LE.0.0) THEN P2 = 0.D0 P1 = 1.D0 END IF Y(IPIX) = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2 XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) ELSE IORD1 = IORD1 + 1 IF (IORD1.GT.NPIXA2) RETURN IORD2 = IORD2 + 1 IF (IORD2.GT.NPIXA2) THEN W0 = 1.E20 ELSE C W0 = WI(IORD2) + DEL W0 = WAVE1(IORD2) END IF C W1 = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL W1 = WAVE2(IORD1) IPIX1 = NINT((WL-WI(IORD1))/YSTEP) + 1 Y(IPIX) = X(IPIX1,IORD1) XMAX = AMAX1(Y(IPIX),XMAX) XMIN = AMIN1(Y(IPIX),XMIN) END IF 20 CONTINUE RETURN END