C @(#)sperebin.for 17.1.1.1 (ES0-DMD) 01/25/02 17:56:15 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 SPEREB C C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 17:23 - 6 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.IDENTIFICATION C C program SPEREBIN version 2.0 840314 C J.D.Ponz ESO - Graching C C.KEYWORDS C C SPECTRA, REBIN C C.PURPOSE C C CONVERT FROM PIXEL INTO WAVELENGTH SAMPLING SPACE C C.ALGORITHM C C EACH LINE IN THE INPUT FRAME IS REBINNED IN LINEAR STEP IN WAVELENGTH C ACCORDING TO THE DISPERSION COEFFS. C C.INPUT/OUTPUT C C REBIN/ECHELLE INPUT OUTPUT STEP [METHOD] C C METHOD CAN BE EITHER "NONLINEAR" - NON LINEAR REBINNING USING C THE COEFFS IN THE TABLE C- IMPLICIT NONE C INTEGER MADRID,KUN,KNUL INTEGER IA1,IAV,II,II1 INTEGER*8 IPNTRA,IPNTRB INTEGER ISTAT INTEGER JJ,NAXISA,NAXISB,NPP,IMIN,IMOUT INTEGER NPIXA(2),NPIXB(2),NCOE(2) INTEGER IMF,TID C REAL WSTART,WSTEP REAL STEPA(2),STARTA(2),STEPB(2),STARTB(2) REAL CUT(4),XMIN,XMAX C DOUBLE PRECISION A1(10),DSTEP(2),DSTART(2) C CHARACTER*80 OUTFR,INFR CHARACTER*8 REDSTA CHARACTER HISTOR*80 CHARACTER IDENTA*72,CUNITA*64,CUNITB*64 CHARACTER*60 TABLE CHARACTER*80 REFER,TABLEF C INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C DATA CUNITB/'FLUX WAVELENGTH'/ C C ... INITIALIZE SYSTEM C CALL STSPRO('SPEREBIN') C C ... READ PARAMETERS C CALL STKRDC('P1',1,1,80,IAV,INFR,KUN,KNUL,ISTAT) CALL STKRDC('P2',1,1,80,IAV,OUTFR,KUN,KNUL,ISTAT) CALL STKRDC('P3',1,1,80,IAV,REFER,KUN,KNUL,ISTAT) II = INDEX('.0123456789',REFER(1:1)) IF (II.EQ.0) THEN CALL STFOPN(REFER,D_R4_FORMAT,0,F_IMA_TYPE,IMF,ISTAT) CALL STDRDD(IMF,'START',1,1,IAV,DSTART,KUN,KNUL,ISTAT) CALL STDRDD(IMF,'STEP',1,1,IAV,DSTEP,KUN,KNUL,ISTAT) CALL STDRDI(IMF,'NPIX',1,1,IAV,NPP,KUn,KNUL,ISTAT) WSTART = DSTART(1) WSTEP = DSTEP(1) ELSE READ (REFER,9000) WSTEP WSTART = -1000. END IF C C CALL STKRDR('INPUTR',1,1,IAV,WSTEP,KUN,KNUL,STAT) CALL STKRDC('MID$CMND',1,1,72,IAV,HISTOR,KUN,KNUL,ISTAT) CALL STKRDC('IN_A',1,1,60,IAV,TABLE,KUN,KNUL,ISTAT) II1 = INDEX(TABLE,' ') - 1 TABLEF = TABLE(1:II1)//'.TBL' CALL TBTOPN(TABLEF,F_I_MODE,TID,ISTAT) CALL STDRDI(TID,'COEFS',1,2,IAV,NCOE,KUN,KNUL,ISTAT) CALL STDRDD(TID,'COEFD',1,10,IAV,A1,KUN,KNUL,ISTAT) CALL STDRDI(TID,'COEFI',1,1,IAV,IA1,KUN,KNUL,ISTAT) DO 10 JJ = IA1 + 1,10 A1(JJ) = 0.D0 10 CONTINUE C C ... MAP INPUT FRAME C CALL STIGET(INFR,10,0,1,2,NAXISA,NPIXA,DSTART,DSTEP, + IDENTA,CUNITA,IPNTRA,IMIN,ISTAT) STARTA(1) = DSTART(1) STARTA(2) = DSTART(2) STEPA(1) = DSTEP(1) STEPA(2) = DSTEP(2) IF (NAXISA.EQ.1 .OR. NPIXA(2).EQ.1) THEN NPIXA(2) = 1 NAXISA = 1 ELSE NAXISA = 2 END IF C C ... MAP OUTPUT FRAME C IF (WSTART.LE.0.0) THEN CALL WRANGE(A1,IA1,WSTEP,NPIXA(1),STARTB(1),NPIXB(1)) ELSE STARTB(1) = WSTART NPIXB(1) = NPP END IF NAXISB = NAXISA NPIXB(2) = NPIXA(2) STARTB(2) = STARTA(2) STEPB(1) = WSTEP STEPB(2) = STEPA(2) REDSTA = 'EXTRACT' DSTART(1) = STARTB(1) DSTART(2) = STARTB(2) DSTEP(1) = STEPB(1) DSTEP(2) = STEPB(2) CALL STIPUT(OUTFR,10,1,1,NAXISB,NPIXB,DSTART,DSTEP, + IDENTA,CUNITB,IPNTRB,IMOUT,ISTAT) C C ... NON LINEAR REBINNIG C CALL REBINB(MADRID(IPNTRA),NAXISA,NPIXA(1),NPIXA(2), + STARTA,STEPA,A1,IA1,MADRID(IPNTRB),NPIXB(1), + NPIXB(2),STARTB,STEPB) C C ... WRITE PROCESS DESCRIPTORS C CALL MNMX(MADRID(IPNTRB),NPIXB(1)*NPIXB(2),XMIN,XMAX) CUT(1) = XMIN CUT(2) = XMAX CUT(3) = XMIN CUT(4) = XMAX CALL STDCOP(IMIN,IMOUT,3,' ',ISTAT) CALL STDWRC(IMOUT,'REDSTAT',1,REDSTA,1,8,KUN,ISTAT) CALL STDWRC(IMOUT,'HISTORY',1,HISTOR,-1,72,KUN,ISTAT) CALL STDWRR(IMOUT,'LHCUTS',CUT,1,4,KUN,ISTAT) DSTART(1) = STARTB(1) DSTART(2) = STARTB(2) CALL STDWRD(IMOUT,'START',DSTART,1,2,KUN,ISTAT) C WRITE (*,*) STARTB(1) CALL STSEPI STOP 9000 FORMAT (F8.4) END DOUBLE PRECISION FUNCTION POLEV(IA,A,XX) C C C IMPLICIT NONE INTEGER IA DOUBLE PRECISION A(IA),XX,YY YY = A(1) + XX* (A(2)+XX* (A(3)+XX* (A(4)+XX* (A(5)+XX* (A(6)+ + XX*A(7)))))) POLEV = YY RETURN END SUBROUTINE MNMX(X,N,XMIN,XMAX) C C FIND MIN MAX VALS C IMPLICIT NONE INTEGER N,I REAL XMIN,XMAX REAL X(N) XMIN = X(1) XMAX = XMIN DO 10 I = 2,N XMIN = AMIN1(XMIN,X(I)) XMAX = AMAX1(XMAX,X(I)) 10 CONTINUE RETURN END SUBROUTINE WRANGE(A,IA,WSTEP,NPIXA,WSTART,NPIXB) C C GET STARTING POSITIONS AND NO. OF PIXELS C IMPLICIT NONE INTEGER IA,NPIXA,NPIXB REAL WSTEP DOUBLE PRECISION WW,A(1),POLEV REAL WSTART C WW = POLEV(IA,A,1.D0) WSTART = WW WW = NPIXA WW = POLEV(IA,A,WW) NPIXB = (WW-WSTART)/WSTEP + 1 RETURN END SUBROUTINE LREBI2(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE LREBI2 VERSION 1 22 MAR 1983 C J.D.PONZ C C.KEYWORDS: C SPECTROSCOPY, CASPEC, REBINNING C C.PURPOSE: C DATA REBINNING C CONVERSION FROM SAMPLES INTO WAVELENGTHS C C.ALGORITHM: C LINEAR INTERPOLATION. THE TRANSFORMATION SAMPLE NUMBER INTO C WAVELENGTH IS DEFINED BY C W = A0 + A1*S + A2*S**2 + ... C C.INPUT/OUTPUT C INPUT ARGUMENTS C C X(NX) REAL*4 INPUT ARRAY C NX INTG*4 NO OF SAMPLES C XSTR REAL*4 STARTING VALUE OF INPUT C XSTP REAL*4 INPUT STEP C AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION C NC REAL*4 NO. OF COEFFICIENTS C NY INTG*4 NO. OF OUTPUT VALUES C YSTR REAL*4 STARTING VALUE OF THE INPUT C YSTP REAL*4 STEP OF THE OUTPUT C C OUTPUT ARGUMENTS C C Y(NY) REAL*4 OUTPUT ARRAY C C C IMPLICIT NONE INTEGER NX,NY,NC,M1,N,MPT,MST,K,KZ6904 REAL XSTR,XSTP,YSTR,YSTP,DELT REAL B2,B3,B4,B5,B6,STARTL,START REAL X(NX),Y(NY) DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,A6,FN,FL,WPOLY DOUBLE PRECISION FACTOR,DPOLY,CNORM DOUBLE PRECISION PN,PNP1,TM,TMP1,C,COV,VAL C WPOLY(FL) = A0 + FL* (A1+FL* (A2+FL* (A3+FL* (A4+FL* (A5+FL* + A6))))) DPOLY(FL) = A1 + FL* (B2+FL* (B3+FL* (B4+FL* (B5+FL*B6)))) C C ... ASSIGN COEFFS C A0 = AC(1) A1 = AC(2) A2 = AC(3) A3 = AC(4) A4 = AC(5) A5 = AC(6) A6 = AC(7) B2 = 2.D0*A2 B3 = 3.D0*A3 B4 = 4.D0*A4 B5 = 5.D0*A5 B6 = 6.D0*A6 C CNORM = 10.D0/YSTP CNORM = 1.D0/YSTP C C ... CHECK FOR OUTPUT START C START = YSTR - 0.5*YSTP STARTL = START M1 = 1 FN = 0.5D0 C PN = 10.D0*WPOLY(FN) PN = WPOLY(FN) DELT = PN - STARTL IF (DELT.GT.0.) THEN M1 = DELT/YSTP + 1. TMP1 = STARTL + M1*YSTP TM = TMP1 - YSTP ELSE TMP1 = STARTL + YSTP TM = STARTL END IF C C ... CHECK FOR INPUT START C DO 10 N = 1,NX FN = FN + 1.D0 PNP1 = WPOLY(FN) C PNP1 = 10.D0*WPOLY(FN) IF (PNP1.GT.TM) GO TO 20 PN = PNP1 10 CONTINUE N = NX + 1 GO TO 30 20 C = X(N) 30 VAL = 0.D0 COV = PNP1 - PN C C ... REBINNING LOOP C C MMSP = NSP MPT = 1 MST = 1 K = 0 99 KZ6904 = 1 IF ( .NOT. (MPT.LE.NY)) GO TO 80 K = K + 1 IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60 C ? 40 IF (PNP1.EQ.PN) GO TO 90 VAL = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV IF (TMP1.LT.PNP1) GO TO 50 N = N + 1 IF (N.GT.NX) GO TO 60 PN = PNP1 FN = FN + 1.D0 C PNP1 = 10.D0*WPOLY(FN) PNP1 = WPOLY(FN) C = X(N) COV = PNP1 - PN GO TO 40 50 TM = TMP1 TMP1 = TMP1 + YSTP 60 FACTOR = CNORM*DPOLY(FN) Y(MPT) = VAL*FACTOR MPT = MPT + 1 VAL = 0. GOTO 99 80 CONTINUE 90 CONTINUE DO 100 K = MPT + 1,NY Y(K) = 0. 100 CONTINUE RETURN END SUBROUTINE REBINB(X,NAXISA,NPIXA1,NPIXA2,STARTA,STEPA, + A1,IA1,Y,NPIXB1,NPIXB2,STARTB,STEPB) C C intermediate routine to rebin each line C IMPLICIT NONE INTEGER NAXISA,NPIXA1,NPIXA2,IA1,NPIXB1,NPIXB2,I REAL STARTA(2),STEPA(2) REAL STARTB(2),STEPB(2) DOUBLE PRECISION A1(1) REAL X(NPIXA1,NPIXA2),Y(NPIXB1,NPIXB2) DO 10 I = 1,NPIXA2 CALL LREBI2(X(1,I),NPIXA1,STARTA,STEPA,A1,IA1, + Y(1,I),NPIXB1,STARTB,STEPB) 10 CONTINUE RETURN END