C @(#)intfr.for 17.1.1.1 (ES0-DMD) 01/25/02 17:11:36 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 SUBROUTINE INTFRI(VAL,NPL,IMAGE,SX,XS,NCUR,NDEG) C+++ C--- IMPLICIT NONE INTEGER NPL REAL VAL(NPL) REAL IMAGE(4) REAL SX REAL XS INTEGER NCUR INTEGER NDEG C INTEGER MAXC INTEGER MAXT PARAMETER (MAXC=100) PARAMETER (MAXT=100000) INTEGER IBAR INTEGER BNMOD INTEGER IAC INTEGER ISTAT INTEGER I INTEGER INNER INTEGER KA INTEGER KUN INTEGER LINET2, LINET3 INTEGER MADRID(1) INTEGER NPL1,NPL2,NPL3 INTEGER NP1, NP2, NP3, NP4 INTEGER NPIX INTEGER NPT REAL AEQW REAL ALICO3, ACONT3 REAL ALICO4, ACONT4 REAL ALINE, ALICO, ACONT REAL ALINT, AITKEN REAL FLXFRC REAL OUTP(10) REAL RC3, XC3, DX3 REAL RC4, XC4, DX4 REAL SX3, SX4 REAL XARG REAL XCUR, YCUR REAL XRANG REAL XPC(MAXC),YPC(MAXC) REAL X(MAXT),Y(MAXT) REAL WORK(MAXT) CHARACTER*80 TEXT,HEADER CHARACTER BIN*4 LOGICAL FIRST INCLUDE 'MID_INCLUDE:PLTDEC.INC/NOLIST' COMMON /VMR/MADRID C DATA IBAR/1/ DATA LINET2/2/ DATA LINET3/3/ C 9010 FORMAT(1X,5G13.6) C C *** start of executable code ********************************************** C Now do the looping for cursers and interpolations C NCUR positions are asked for NPL1 = INT(MIN(IMAGE(1),IMAGE(2))) NPL2 = INT(MAX(IMAGE(1),IMAGE(2))) NPL3 = 1 NPT = ABS(NPL1-NPL2) + 1 C C *** Now do the looping for cursers and integration C Actions: NCUR cursors are asked for C A space bar terminates this loop FIRST = .TRUE. 30 CONTINUE DO 50 I = 1,NCUR CALL PTGCUR(XCUR,YCUR,KA,ISTAT) IF (KA.EQ.32) THEN ! exit on space bar GO TO 120 ELSE ! correct input; store data point XPC(I) = XCUR YPC(I) = YCUR CALL PTDATA(5,0,0,XPC(I),YPC(I),0.0,1) END IF 50 CONTINUE CALL SORT2(XPC,YPC,NCUR) C C *** Determine the range effective for change (= the EXTREME cursors) NP1 = NINT((XPC(1)-XS)/SX) + 1 ! pixel number corr. to cursor NP2 = NINT((XPC(NCUR)-XS)/SX) + 1 IF (NP1.LT.NPL1) THEN NP1 = NPL1 ENDIF C C *** Do the interpolation here and plot the line IF (NDEG.LT.2) THEN DO 70 I = NP1,NP2 XARG = XS + FLOAT(I-1)*SX WORK(I) = ALINT(XPC,YPC,NCUR,XARG,1,NCUR) 70 CONTINUE ELSE DO 80 I = NP1,NP2 XARG = XS + FLOAT(I-1)*SX WORK(I) = AITKEN(XPC,YPC,NCUR,XARG,1,NCUR,NDEG) 80 CONTINUE END IF C C *** get the bin option, plot and write the results CALL PTKRDC('BINMOD',4,IAC,BIN) IF (BIN(1:2).EQ.'ON') THEN BNMOD = 1 ELSE BNMOD = 0 ENDIF NPIX = NP2 - NP1 + 1 DO 90 I = NP1,NP2 X(I) = XS + FLOAT(I-1)*SX Y(I) = WORK(I) 90 CONTINUE CALL PTDATA(0,LINET3,BNMOD,X(NP1),Y(NP1),0.0,NPIX) C C *** Do the integration C FIRST compute the flux density in the first bin. C The bin limits are the input cursor position and the middle between C the first and second image pixel. We compute the interpolated values C at these points and make the integral in the first bin. INNER = NCUR/2 ! first inner pixel XC3 = XPC(INNER) ! real cursor position RC3 = (XC3-XS)/SX+1.0 ! real cursor pixel NP3 = NINT((XC3-XS)/SX)+1 ! nearest image pixel DX3 = FLOAT(NP3)-RC3 SX3 = (DX3+0.5)*SX ACONT3 = SX3*WORK(NP3) ALICO3 = SX3*VAL(NP3) C C *** SECOND compute the flux density in the last bin. C The bin limits are the last cursor position and the middle between C the last and one but last image pixel. We compute the interpolated C values at these points and make the integral in the last bin. C and then do the integration. XC4 = XPC(INNER+1) RC4 = (XC4-XS)/SX+1.0 NP4 = NINT((XC4-XS)/SX)+1 DX4 = RC4-FLOAT(NP4) SX4 = (DX4+0.5)*SX ACONT4 = SX4*WORK(NP4) ALICO4 = SX4*VAL(NP4) C C *** do the pixel summation for continuum and line ACONT = 0. ALICO = 0. DO 110 I = NP3+1,NP4-1 ACONT = WORK(I)*SX + ACONT ALICO = VAL(I)*SX + ALICO 110 CONTINUE ACONT = ACONT + ACONT3 + ACONT4 ALICO = ALICO + ALICO3 + ALICO4 XRANG = XPC(INNER+1)-XPC(INNER) ALINE = ALICO - ACONT AEQW = -1.0*(ALINE/ACONT)*XRANG FLXFRC = ALINE/ACONT C *** C IF (FIRST) THEN HEADER = ' X_start (pix/world) '// 2 ' X_end (pix/world) Pixel sep. ' CALL STTPUT(HEADER,ISTAT) HEADER = ' Line+Cont. Continuum '// 2 'Line Line/Cont Equiv. w. ' CALL STTPUT(HEADER,ISTAT) HEADER = ' ---------------------------------------'// 2 '-------------------------' CALL STTPUT(HEADER,ISTAT) FIRST = .FALSE. END IF WRITE(TEXT,9010) XC3,RC3,XC4,RC4,SX CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,9010) ALICO,ACONT,ALINE,FLXFRC,AEQW CALL STTPUT(TEXT,ISTAT) C *** C OUTP(1) = RC3 ! first NUMBER OUTP(2) = RC4 ! last pixel OUTP(3) = XC3 ! first cursor OUTP(4) = XC4 ! last cursor OUTP(5) = SX ! step in x OUTP(6) = ALICO ! line + cont OUTP(7) = ACONT ! continuum OUTP(8) = ALINE ! line OUTP(9) = FLXFRC ! flux fraction OUTP(10) = AEQW ! equivalent width CALL STKWRR('OUTPUTR',OUTP,1,10,KUN,ISTAT) GO TO 30 120 CONTINUE RETURN END SUBROUTINE INTFRB(VAL,NPL,IMAGE,SX,XS,NCUR,NDEG, 2 DELTAX,XPC) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.AUTHOR: F. Murtagh, ST-ECF C.VERSION: 860321 FM Creation C.VERSION: 880607 RHW Restructured and recoded for portable MIDAS C.VERSION: 900430 RHW Make code similar to INTFRI C.VERSION: 900430 RHW Store more output in the keyword OUTPUTR C.VERSION: 900629 RHW Integration over subpixel included C.VERSION: 910115 RHW IMPLICIT NONE added C---------------------------------------------------------------------------- IMPLICIT NONE INTEGER MAXC PARAMETER (MAXC=100) C INTEGER NPL REAL VAL(NPL) REAL IMAGE(4) REAL SX REAL XS INTEGER NCUR INTEGER NDEG REAL DELTAX REAL XPC(MAXC) C INTEGER MAXT PARAMETER (MAXT=100000) INTEGER ISTAT INTEGER KUN INTEGER I, I2, INNER INTEGER NP1, NP2, NP3, NP4 INTEGER NPL1, NPL2, NPL3, NPT REAL AITKEN, ALINE, ALINT REAL ACONT, ALICO REAL ACONT3, ACONT4, ALICO3, ALICO4 REAL AEQW REAL DX3, DX4 REAL FLXFRC REAL OUTP(10) REAL RC3, RC4 REAL SX3, SX4 REAL WORK(MAXT) REAL XARG REAL Y REAL XC3, XC4 REAL XRANG REAL YPC(MAXC) C C *** start of executable code ********************************************** C Now do the looping for cursers and interpolations. C NCUR positions are asked for. NPL1 = INT(MIN(IMAGE(1),IMAGE(2))) NPL2 = INT(MAX(IMAGE(1),IMAGE(2))) NPL3 = 1 NPT = ABS(NPL1-NPL2) + 1 C C *** loop through cursor position and determine nearest pixels DO 10 I = 1, NCUR NP1 = NINT((XPC(I)-DELTAX-XS)/SX) + 1 NP2 = NINT((XPC(I)+DELTAX-XS)/SX) + 1 IF (NP1.LT.NPL1 .OR. NP1.GT.NPL2 .OR. 2 NP2.LT.NPL1 .OR. NP2.GT.NPL2) THEN CALL STTPUT('*** FATAL: Cursor position(s) outside'// 2 ' frame; try again', ISTAT) CALL STSEPI ENDIF Y = 0.0 DO 20 I2 = NP1, NP2 Y = Y + VAL(I2) 20 CONTINUE Y = Y/FLOAT(NP2-NP1+1) YPC(I) = Y 10 CONTINUE CALL SORT2(XPC,YPC,NCUR) C C *** Determine effective range of cursor positions NP1 = NINT((XPC(1)-XS)/SX) + 1 ! first cursor position NP2 = NINT((XPC(NCUR)-XS)/SX) + 1 ! last cursor position C C *** Do the interpolation here IF (NDEG.LT.2) THEN !compute the continuum DO 30 I = NP1,NP2 XARG = XS + FLOAT(I-1)*SX WORK(I) = ALINT(XPC,YPC,NCUR,XARG,1,NCUR) 30 CONTINUE ELSE DO 40 I = NP1,NP2 XARG = XS + FLOAT(I-1)*SX WORK(I) = AITKEN(XPC,YPC,NCUR,XARG,1,NCUR,NDEG) 40 CONTINUE ENDIF C C *** Do the integration C FIRST compute the flux density in the first bin. C The bin limits are the input cursor position and the middle between C the first and second image pixel. We compute the interpolated values C at these points and make the integral in the first bin. INNER = NCUR/2 ! first inner pixel XC3 = XPC(INNER) ! real cursor position RC3 = (XC3-XS)/SX+1.0 ! real cursor pixel NP3 = NINT((XC3-XS)/SX)+1 ! nearest image pixel DX3 = FLOAT(NP3)-RC3 SX3 = (DX3+0.5)*SX ACONT3 = SX3*WORK(NP3) ALICO3 = SX3*VAL(NP3) C C *** SECOND compute the flux density in the last bin. C The bin limits are the last cursor position and the middle between C the last and one but last image pixel. We compute the interpolated C values at these points and make the integral in the last bin. C and then do the integration. XC4 = XPC(INNER+1) RC4 = (XC4-XS)/SX+1.0 NP4 = NINT((XC4-XS)/SX)+1 DX4 = RC4-FLOAT(NP4) SX4 = (DX4+0.5)*SX ACONT4 = SX4*WORK(NP4) ALICO4 = SX4*VAL(NP4) C C *** do the pixel summation for continuum and line ACONT = 0. ALICO = 0. DO 50 I = NP3+1,NP4-1 ACONT = WORK(I)*SX + ACONT ALICO = VAL(I)*SX + ALICO 50 CONTINUE ACONT = ACONT + ACONT3 + ACONT4 ALICO = ALICO + ALICO3 + ALICO4 XRANG = XPC(INNER+1)-XPC(INNER) ALINE = ALICO - ACONT AEQW = -1.0*(ALINE/ACONT)*XRANG FLXFRC = ALINE/ACONT C *** C OUTP(1) = RC3 ! first NUMBER OUTP(2) = RC4 ! last pixel OUTP(3) = XC3 ! first cursor OUTP(4) = XC4 ! last cursor OUTP(5) = SX ! step in x OUTP(6) = ALICO ! line + cont OUTP(7) = ACONT ! continuum OUTP(8) = ALINE ! line OUTP(9) = FLXFRC ! flux fraction OUTP(10) = AEQW ! equivalent width CALL STKWRR('OUTPUTR',OUTP,1,10,KUN,ISTAT) C RETURN END