SUBROUTINE ZPPF * * Module Number: 13.6 * * Module Name: ZPPF * * Keyphrase: * ---------- * Compute HRS paired pulse coefficients * * Description: * ------------ * The paired pulse equation is given by: * * y = (1-exp(-tx)/t * where * x is the true count rate * y is the observed count rate * t is a time coef. * * This routine computes the value of t using the * Levenberg-Marquardt non-linear least squares method * from NUMERCIAL RECIPES. * * Fortran Name: ZPPF * * Keywords of accessed files and tables: * -------------------------------------- * OBSERVED input Observed count rate file * EXPECTED input Expected count rate file * TABLE output table of fit results * PPTAB output table of paired pulse coef. for RSDP * TAU input initial guess for time coef. * NITER input Number of least squares iterations * * * Subroutines Called: * ------------------- * CDBS: * ZPPF1 * SDAS: * uclgs* , umsput * uttinn, utppti, utcdef, utrpt*, utcpt*, uthad*, uttclo, utccre * uttopn, utpgti, utcfnd, utrgt*, utcgt*, uthgt*, uttclo * uimotp, uimxtp, uimctp, uimopn, uimgid, uimclo, uhdgs* uimclo * * History: * -------- * version date Author Description * 1 4/20/87 D. Lindler Designed and coded * 2 Aug 1988 D. Lindler New paired pulse formula * switched to non-linear method *------------------------------------------------------------------------- C INCLUDE FILE FOR THE IRAF77 FORTRAN INTERFACE TO THE IRAF VOS C C C FILE I/O ACCESS MODES C INTEGER RDONLY PARAMETER (RDONLY = 1) INTEGER RDWRIT PARAMETER (RDWRIT = 2) INTEGER WRONLY PARAMETER (WRONLY = 3) INTEGER APPEND PARAMETER (APPEND = 4) C C CODES FOR DATA TYPES C INTEGER TYBOOL PARAMETER (TYBOOL = 1) INTEGER TYCHAR PARAMETER (TYCHAR = 2) INTEGER TYINT PARAMETER (TYINT = 4) INTEGER TYREAL PARAMETER (TYREAL = 6) INTEGER TYDOUB PARAMETER (TYDOUB = 7) C C UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87 C INTEGER STDOUT PARAMETER (STDOUT = 1) INTEGER STDERR PARAMETER (STDERR = 2) C C UHDAS HEADER PARM TYPES -- CB, DAO, 5-SEP-87 C INTEGER GENHDR PARAMETER (GENHDR = 0) INTEGER IMSPEC PARAMETER (IMSPEC = 1) C C THIS SECTION IS FOR PARAMETERS RELEVANT TO TABLE I/O. C C THESE MAY BE SET BY UTPPTI AND/OR READ BY UTPGTI: C C LENGTH OF ROW (UNIT = SIZE OF REAL) INTEGER TBRLEN PARAMETER (TBRLEN = 1) C INCREASE ROW LENGTH INTEGER TBIRLN PARAMETER (TBIRLN = 2) C NUMBER OF ROWS TO ALLOCATE INTEGER TBALLR PARAMETER (TBALLR = 3) C INCREASE ALLOC NUM OF ROWS INTEGER TBIALR PARAMETER (TBIALR = 4) C WHICH TYPE OF TABLE? (ROW OR COLUMN) INTEGER TBWTYP PARAMETER (TBWTYP = 5) C MAXIMUM NUMBER OF USER PARAMETERS INTEGER TBMXPR PARAMETER (TBMXPR = 6) C MAXIMUM NUMBER OF COLUMNS INTEGER TBMXCL PARAMETER (TBMXCL = 7) C TYPE = ROW-ORDERED TABLE INTEGER TBTYPR PARAMETER (TBTYPR = 11) C TYPE = COLUMN-ORDERED TABLE INTEGER TBTYPC PARAMETER (TBTYPC = 12) C C THESE MAY BE READ BY UTPGTI BUT MAY NOT BE SET: C C NUMBER OF ROWS WRITTEN TO INTEGER TBNROW PARAMETER (TBNROW = 21) C C END IRAF77.INC C C ERROR PROCESSING PARAMETERS C CHARACTER*130 CONTXT INTEGER STATUS,ISTAT,ISTATS(20) C C OUTPUT PAIRED PULSE COEF. TABLE C CHARACTER*64 PPTAB INTEGER IDOUT,COLIDS(9),CTYPE(9),NROWS CHARACTER*12 COLNAM(9),CFORM(9),CUNITS(9) C C OUTPUT RESULT TABLE C CHARACTER*64 TABLE CHARACTER*15 CNAME2(4) INTEGER CTYPE2(4) C C INPUT FILE I/O C CHARACTER*130 OBSER,EXPECT,NAME,NAME2 INTEGER IDTEMP,IDTMP2,IDIN,IDIN2,NAXIS,DTYPE,DIMEN(8) INTEGER NS,NPOINT C C INPUT KEYWORD PARAMETERS C DOUBLE PRECISION TAU,THRESH INTEGER IDET,DET1,DET2,NITER C C FIT PARAMETERS C DOUBLE PRECISION OBS(50000),EX(50000) C C OTHER LOCAL VARIABLES C INTEGER I C C DATA DECLARATIONS C DATA COLNAM/'INSTRUMENT','DETECTOR','TAU1','EPSILON','THRESHOLD' * ,'Q0','Q1','F','ITERATIONS'/ DATA CTYPE/-3,TYINT,6*TYDOUB,TYINT/ DATA CNAME2/'EXPECTED','OBSERVED','OBS_FIT','RESIDUAL'/ DATA CTYPE2/4*TYREAL/ C-------------------------------------------------------------------------- C C READ CL PARAMETERS C CALL UCLGST('observed',OBSER,ISTATS(1)) CALL UCLGST('expected',EXPECT,ISTATS(2)) CALL UCLGST('table',TABLE,ISTATS(3)) CALL UCLGST('pptab',PPTAB,ISTATS(4)) CALL UCLGSD('tau',TAU,ISTATS(5)) CALL UCLGSI('niter',NITER,ISTATS(6)) CALL UCLGSD('threshold',THRESH,ISTATS(7)) CALL UCLGSI('detector',IDET,ISTATS(8)) DO 10 I=1,8 IF(ISTATS(I).NE.0)THEN CONTXT='Error reading CL parameter' GO TO 999 ENDIF 10 CONTINUE C C READ INPUT DATA --------------------------------------------------------- C NPOINT=0 C --->TOTAL NUMBER OF DATA POINTS READ C C OPEN TEMPLATES C CALL UIMOTP(OBSER,IDTEMP,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening observed rate file template '// * OBSER GO TO 999 ENDIF CALL UIMOTP(EXPECT,IDTMP2,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening expected rate file template '// * EXPECT GO TO 999 ENDIF C C GET NEXT FILE NAME C 20 CALL UIMXTP(IDTEMP,NAME,ISTAT) IF(ISTAT.LT.0)GO TO 100 IF(ISTAT.NE.0)THEN CONTXT='Error getting filename from template '//OBSER GO TO 999 ENDIF C C GET NEXT EXPECTED COUNT RATE FILE NAME C CALL UIMXTP(IDTMP2,NAME2,ISTAT) IF(ISTAT.LT.0)THEN CONTXT='Insufficent expected rate files specified' GO TO 999 ENDIF IF(ISTAT.NE.0)THEN CONTXT='Error getting filename from template '//EXPECT GO TO 999 ENDIF C C OPEN INPUT FILES C CALL UIMOPN(NAME,RDONLY,IDIN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening input file '//NAME GO TO 999 ENDIF CALL UIMOPN(NAME2,RDONLY,IDIN2,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening input file '//NAME2 GO TO 999 ENDIF C C READ IMAGE INFO C CALL UIMGID(IDIN,DTYPE,NAXIS,DIMEN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//NAME GO TO 999 ENDIF C C CHECK FOR VALID DATA C IF(NAXIS.NE.1)THEN CONTXT='Input data must be one dimensional' GO TO 999 ENDIF NS=DIMEN(1) C C CHECK TO SEE IF WE STILL HAVE ROOM C IF((NS+NPOINT).GT.50000)THEN CONTXT='Total number of input data points can not '// * 'exceed 50000' GO TO 999 ENDIF C C READ EXPECTED COUNT RATE FILE FILE INFO C CALL UIMGID(IDIN2,DTYPE,NAXIS,DIMEN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading file '//NAME2 GO TO 999 ENDIF C C CHECK FOR VALID DATA C IF((NAXIS.NE.1).OR.(DIMEN(1).NE.NS))THEN CONTXT='Expected count rate file '// * 'must be same size as observed' GO TO 999 ENDIF C C READ DATA C CALL UIGL1D(IDIN,OBS(NPOINT),ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//NAME GO TO 999 ENDIF CALL UIGL1D(IDIN2,EX(NPOINT),ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//NAME2 GO TO 999 ENDIF NPOINT=NPOINT+NS C C CLOSE IMAGES C CALL UIMCLO(IDIN,ISTAT) CALL UIMCLO(IDIN2,ISTAT) C C GO GET NEXT IMAGE C GO TO 20 100 CONTINUE C C DONE READING INPUT DATA ------------------------------------------------ IF(NPOINT.EQ.0)THEN CONTXT='No input files found' GO TO 999 ENDIF C C C COMPUTE PAIRED PULSE COEF. AND FITS C CALL ZPPF1(OBS,EX,NPOINT,NITER,TAU,STATUS) IF(STATUS.NE.0)THEN CONTXT='No output files generated' GO TO 999 ENDIF C C WRITE FIT RESULT TABLE --------------------------------------------------- C C C OPEN OUTPUT TABLE C CALL UTTINN(TABLE,IDOUT,ISTATS(1)) CALL UTPPTI(IDOUT,TBWTYP,TBTYPC,ISTATS(2)) CALL UTPPTI(IDOUT,TBALLR,NPOINT,ISTATS(3)) CALL UTCDEF(IDOUT,CNAME2,CUNITS,CFORM,CTYPE2,4,COLIDS,ISTATS(4)) CALL UTTCRE(IDOUT,ISTATS(5)) DO 200 I=1,5 IF(ISTATS(I).NE.0)THEN CONTXT='Error creating output table '//TABLE GO TO 999 ENDIF 200 CONTINUE C C COPY RESULTS TO TABLE C CALL UTCPTD(IDOUT,COLIDS(1),1,NPOINT,EX,ISTATS(1)) CALL UTCPTD(IDOUT,COLIDS(2),1,NPOINT,OBS,ISTATS(2)) C COMPUTE OBS_FIT AND SAVE IN EX VECTOR DO 201 I=1,NPOINT 201 EX(I)=(1-DEXP(-TAU*EX(I)))/TAU CALL UTCPTD(IDOUT,COLIDS(3),1,NPOINT,EX,ISTATS(3)) C COMPUTE RESIDUALS AND SAVE IN EX VECTOR DO 202 I=1,NPOINT 202 EX(I)=OBS(I)-EX(I) CALL UTCPTD(IDOUT,COLIDS(4),1,NPOINT,EX,ISTATS(4)) DO 210 I=1,4 IF(ISTATS(I).NE.0)THEN CONTXT='Error writing to output table '//TABLE GO TO 999 ENDIF 210 CONTINUE CALL UTTCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output table '//TABLE GO TO 999 ENDIF C C WRITE PAIRED PULSE COEFFICIENT TABLE ------------------------------------ C C C OPEN OUTPUT TABLE C CALL UTTINN(PPTAB,IDOUT,ISTATS(1)) CALL UTPPTI(IDOUT,TBRLEN,15,ISTATS(2)) CALL UTPPTI(IDOUT,TBMXCL,9,ISTATS(3)) CALL UTCDEF(IDOUT,COLNAM,CUNITS,CFORM,CTYPE,9,COLIDS,ISTATS(4)) CALL UTTCRE(IDOUT,ISTATS(5)) DO 300 I=1,5 IF(ISTATS(I).NE.0)THEN CONTXT='Error creating output table '//PPTAB GO TO 999 ENDIF 300 CONTINUE C C COPY RESULTS TO TABLE, ONE ROW FOR EACH DETECTOR IF IDET=0 C OTHERWISE ONE ROW FOR THE SEPCIFIED IDET C IF(IDET.EQ.0)THEN DET1=1 DET2=2 ELSE DET1=IDET DET2=IDET ENDIF NROWS=0 DO 700 IDET=DET1,DET2 NROWS=NROWS+1 CALL UTRPTT(IDOUT,COLIDS(1),1,NROWS,'HRS',ISTATS(1)) CALL UTRPTI(IDOUT,COLIDS(2),1,NROWS,IDET,ISTATS(2)) CALL UTRPTD(IDOUT,COLIDS(3),1,NROWS,TAU,ISTATS(3)) CALL UTRPTD(IDOUT,COLIDS(4),1,NROWS,0.0D0,ISTATS(4)) CALL UTRPTD(IDOUT,COLIDS(5),1,NROWS,THRESH,ISTATS(5)) CALL UTRPTD(IDOUT,COLIDS(6),1,NROWS,0.0D0,ISTATS(6)) CALL UTRPTD(IDOUT,COLIDS(7),1,NROWS,0.0D0,ISTATS(7)) CALL UTRPTD(IDOUT,COLIDS(8),1,NROWS,0.0D0,ISTATS(8)) CALL UTRPTI(IDOUT,COLIDS(9),1,NROWS,1,ISTATS(9)) DO 310 I=1,9 IF(ISTATS(I).NE.0)THEN CONTXT='Error writing to output table'//PPTAB GO TO 999 ENDIF 310 CONTINUE 700 CONTINUE CALL UTTCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output table '//TABLE GO TO 999 ENDIF C C DONE C GO TO 1000 999 CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT) 1000 RETURN END