SUBROUTINE YAFITE * * Module Number: 14.9.2 * * Module Name: YAFITE * * Keyphrase: * ---------- * smooth sensitvity curve * * Description: * ------------ * The routine computes a smooth sensitivity curve by fitting * a least squares cubic spline function of specified number of * nodes to the input curve or by using supplied input node * positions in an input table. * * Fortran Name: zafite.for * * Keywords of accessed files and tables: * -------------------------------------- * table input input raw sens. table * wave input vector of wavelengths to compute inv. sens. for * output output output inverse sens. reference file * nodes input number of spline nodes * append input append fit to input table? * innodes input table of input node positions * outnodes output table of output node/value positions * niter input number of iterations of least squares fit * template input template for output file * * Subroutines Called: * ------------------- * CDBS: * zspfit, zsplin, ztplat, yxptrn * 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 3/11/87 D. Lindler Designed and coded * 2 May 88 D. Lindler New sdas i/o and standards * 2.1 May 94 H. Bushouse Added YBASE,YSPACE to YXPTRN *------------------------------------------------------------------------- C C FILE I/O ACCESS MODES C INTEGER RDONLY PARAMETER (RDONLY = 1) INTEGER RDWRIT PARAMETER (RDWRIT = 2) C C CODES FOR DATA TYPES C 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) INTEGER USRLOG PARAMETER (USRLOG = 4) 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 NUMBER OF COLUMNS DEFINED INTEGER TBNCOL PARAMETER (TBNCOL = 22) C AMOUNT OF ROW USED (UNIT=SIZE OF REAL) INTEGER TBRUSD PARAMETER (TBRUSD = 23) C NUMBER OF USER PARAMETERS INTEGER TBNPAR PARAMETER (TBNPAR = 24) C END IRAF77.INC C C ERROR PROCESSING PARAMETERS C INTEGER STATUS,ISTAT,ISTATS(20) CHARACTER*130 CONTXT C C INPUT CL PARAMETERS C CHARACTER*64 TABLE,INNODE,OUTNOD,TMPLT,WNAME,OUTPUT INTEGER NODES,NITER LOGICAL APPEND C C INPUT FILE I/O C INTEGER IDIN,NAXIS,DTYPE,DIMEN(8),NOUT DOUBLE PRECISION WOUT(3000) C C INPUT TABLE VARIABLES C INTEGER COLIDS(2) CHARACTER*16 COLNAM(2) DOUBLE PRECISION WAVE(3000),VALUES(3000) INTEGER NROWS,NS LOGICAL NULLS(3000) C C OUTPUT NODE TABLE C DOUBLE PRECISION WNODES(200),VNODES(200) CHARACTER*8 CUNITS(2),CFORM(2) INTEGER IDOUT,CTYPE(2) C C FIT ROUTINE PARAMETERS C DOUBLE PRECISION FIT(3000),FIT2(3000) C C OBSERVING MODE PARAMETERS C CHARACTER*5 DET CHARACTER*3 APER,GRAT CHARACTER*6 APERPS CHARACTER*1 POLAR INTEGER PASSDR,XSTEPS INTEGER FCHNL,NCHNLS,OVRSCN,YBASE,YSPACE C C OTHERS C INTEGER I DOUBLE PRECISION DIF C C DATA DECLARATIONS C DATA COLNAM/'WAVELENGTH','VALUE'/ DATA CUNITS/2*' '/ DATA CFORM/2*' '/ DATA CTYPE/2*TYDOUB/ C----------------------------------------------------------------------- C C GET INPUT CL PARAMETERS C CALL UCLGST('table',TABLE,ISTATS(1)) CALL UCLGST('wave',WNAME,ISTATS(2)) CALL UCLGST('output',OUTPUT,ISTATS(3)) CALL UCLGSI('nodes',NODES,ISTATS(4)) CALL UCLGSI('niter',NITER,ISTATS(5)) CALL UCLGSB('append',APPEND,ISTATS(6)) CALL UCLGST('innodes',INNODE,ISTATS(7)) CALL UCLGST('outnodes',OUTNOD,ISTATS(8)) CALL UCLGST('template',TMPLT,ISTATS(9)) DO 10 I=1,9 IF(ISTATS(I).NE.0)THEN CONTXT='Error obtaining CL parameter' GO TO 999 ENDIF 10 CONTINUE C C READ INPUT TABLE OF WAVELENGTH VERSUS RAW SENSITIVITY VALUES---------- C CALL UTTOPN(TABLE,RDONLY,IDIN,ISTAT) C --->open table IF(ISTAT.NE.0)THEN CONTXT='Error opening input table '//TABLE GO TO 999 ENDIF CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT) C --->get number of rows IF(ISTAT.NE.0)THEN CONTXT='Error reading input table '//TABLE GO TO 999 ENDIF IF(NROWS.GT.3000)THEN CONTXT='Max. of 3000 rows allowed in sens. table '//TABLE GO TO 999 ENDIF CALL UTCFND(IDIN,COLNAM,2,COLIDS,ISTAT) C --->locate columns IF(ISTAT.NE.0)THEN CONTXT='Error locating columns in sens. table '//TABLE GO TO 999 ENDIF C C READ TABLE DATA C NS=NROWS CALL UTCGTD(IDIN,COLIDS(1),1,NS,WAVE,NULLS,ISTATS(1)) CALL UTCGTD(IDIN,COLIDS(2),1,NS,VALUES,NULLS,ISTATS(2)) IF((ISTATS(1).NE.0).OR.(ISTATS(2).NE.0))THEN CONTXT='Error reading input table '//TABLE GO TO 999 ENDIF C C GET CONFIGURATION DATA C CALL UTHGTT(IDIN,'DETECTOR',DET,ISTATS(1)) CALL UTHGTT(IDIN,'APER_ID',APER,ISTATS(2)) CALL UTHGTT(IDIN,'FGWA_ID',GRAT,ISTATS(3)) CALL UTHGTT(IDIN,'APER_POS',APERPS,ISTATS(4)) CALL UTHGTT(IDIN,'POLAR_ID',POLAR,ISTATS(5)) CALL UTHGTI(IDIN,'PASS_DIR',PASSDR,ISTATS(6)) DO 90 I=1,10 IF(ISTATS(I).NE.0)THEN CONTXT='ERROR getting configuration data for input table' GO TO 999 ENDIF 90 CONTINUE CALL UTTCLO(IDIN,ISTAT) C C READ INPUT NODE TABLE (IF NOT SUPPLIED) USE NUMBER OF NODES TO COMPUTE ----- C EQUALLY SPACED NODE POSITIONS C IF(INNODE.NE.' ')THEN CALL UTTOPN(INNODE,RDONLY,IDIN,ISTAT) C --->open table IF(ISTAT.NE.0)THEN CONTXT='Error opening input table '//INNODE GO TO 999 ENDIF CALL UTPGTI(IDIN,TBNROW,NODES,ISTAT) C --->get number of rows IF(ISTAT.NE.0)THEN CONTXT='Error reading input table '//INNODE GO TO 999 ENDIF IF((NODES.GT.30).OR.(NODES.LT.2))THEN CONTXT='Error: '//INNODE//' Must have between 2 and '// * '30 node positions' GO TO 999 ENDIF CALL UTCFND(IDIN,COLNAM,1,COLIDS,ISTAT) C --->locate column IF(ISTAT.NE.0)THEN CONTXT='Error locating wavelength column in sens. table ' * //INNODE GO TO 999 ENDIF C C READ TABLE DATA C CALL UTCGTD(IDIN,COLIDS(1),1,NODES,WNODES,NULLS,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input table '//INNODE GO TO 999 ENDIF CALL UTTCLO(IDIN,ISTAT) CALL UMSPUT('Input node positions taken from table '// * INNODE,STDOUT,0,ISTAT) ELSE C C USE EQUALLY SPACED NODES (UNLESS NODES EQUAL 0, WHERE NO FIT IS DONE) C IF(NODES.GT.1)THEN DIF=WAVE(NS)-WAVE(1) DIF=DIF/(NODES-1) DO 20 I=1,NODES WNODES(I)=WAVE(1)+DIF*(I-1) 20 CONTINUE ENDIF ENDIF C C READ WAVELENGTH VECTOR TO USE FOR OUTPUT SENSITIVITIES ---------------------- C CALL UIMOPN(WNAME,RDONLY,IDIN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening wave file '//WNAME 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 '//WNAME GO TO 999 ENDIF C C CHECK FOR VALID DATA C IF((NAXIS.NE.1).OR.(DIMEN(1).GT.3000))THEN CONTXT='Input data must be a vector (up to 3000 points)' GO TO 999 ENDIF NOUT=DIMEN(1) C C GET XSTEP PATTERN INFORMATION C CALL YXPTRN(IDIN,FCHNL,NCHNLS,XSTEPS,OVRSCN,YBASE,YSPACE, * ISTAT) IF(ISTAT.NE.0)THEN CONTXT='PATTERN INFORMATION MISSING FROM '//WNAME GO TO 999 ENDIF C C READ DATA C CALL UIGL1D(IDIN,WOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//WNAME GO TO 999 ENDIF C C CLOSE IMAGE C CALL UIMCLO(IDIN,ISTAT) C C PERFORM LEAST SQUARES SPLINE FIT C CALL ZSPFIT(WAVE,VALUES,NS,WNODES,VNODES,NODES,NITER, * FIT,STATUS) IF(STATUS.NE.0)THEN CONTXT='ERROR DOING LEAST SQUARES SPLINE FIT' GO TO 999 ENDIF C C COMPUTE SPLINE CURVE AT POINTS IN WAVE C CALL ZSPLIN(WOUT,NOUT,WNODES,VNODES,NODES,FIT2,STATUS) IF(STATUS.NE.0)THEN CONTXT='ERROR COMPUTING SPLINE FUNCTION' GO TO 999 ENDIF C C CONVERT TO INVERSE SENSITIVITY C DO 700 I=1,NOUT IF((WOUT(I).EQ.0.0).OR.(FIT2(I).LE.0))THEN FIT2(I)=0.0D0 ELSE FIT2(I)=1.0D0/FIT2(I) ENDIF 700 CONTINUE C C WRITE OUTPUT TABLE OF NODES ------------------------------------------- C CALL UTTINN(OUTNOD,IDOUT,ISTATS(1)) CALL UTPPTI(IDOUT,TBWTYP,TBTYPC,ISTATS(2)) CALL UTPPTI(IDOUT,TBMXCL,2,ISTATS(3)) CALL UTPPTI(IDOUT,TBALLR,NODES,ISTATS(4)) CALL UTCDEF(IDOUT,COLNAM,CUNITS,CFORM,CTYPE,2,COLIDS,ISTATS(5)) CALL UTTCRE(IDOUT,ISTATS(6)) DO 200 I=1,6 IF(ISTATS(I).NE.0)THEN CONTXT='Error creating output table '//OUTNOD GO TO 999 ENDIF 200 CONTINUE C C COPY RESULTS TO TABLE C CALL UTCPTD(IDOUT,COLIDS(1),1,NODES,WNODES,ISTATS(1)) CALL UTCPTD(IDOUT,COLIDS(2),1,NODES,VNODES,ISTATS(2)) DO 210 I=1,2 IF(ISTATS(I).NE.0)THEN CONTXT='Error writing to output node table' GO TO 999 ENDIF 210 CONTINUE CALL UTTCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output table '//OUTNOD GO TO 999 ENDIF C C WRITE OUTPUT REFERENCE FILE USING INPUT TEMPLATE FILES C CALL ZTPLAT(TMPLT,OUTPUT,FIT2,NOUT,IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Unable to write output inv. sensitivity file' GO TO 999 ENDIF CALL UHDPST(IDOUT,'FGWA_ID',GRAT,ISTATS(1)) CALL UHDPST(IDOUT,'DETECTOR',DET,ISTATS(2)) CALL UHDPST(IDOUT,'APER_ID',APER,ISTATS(3)) CALL UHDPST(IDOUT,'POLAR_ID',POLAR,ISTATS(4)) CALL UHDPST(IDOUT,'APER_POS',APERPS,ISTATS(5)) CALL UHDPSI(IDOUT,'PASS_DIR',PASSDR,ISTATS(6)) CALL UHDPSI(IDOUT,'FCHNL',FCHNL,ISTATS(7)) CALL UHDPSI(IDOUT,'NCHNLS',NCHNLS,ISTATS(8)) CALL UHDPSI(IDOUT,'NXSTEPS',XSTEPS,ISTATS(9)) CALL UHDPSI(IDOUT,'OVERSCAN',OVRSCN,ISTATS(10)) DO 910 I=1,10 IF(ISTATS(I).NE.0)THEN CONTXT='ERROR updating keyword value in '// * 'output inv. sens. file' ENDIF 910 CONTINUE CALL UIMCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing file '//OUTPUT GO TO 999 ENDIF C C ADD FITTED VALUES TO INPUT TABLE IF APPEND IS SPECIFIED----------------------- C IF(APPEND)THEN CALL UTTOPN(TABLE,RDWRIT,IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='ERROR opening table for output '// * TABLE GO TO 999 ENDIF C C CHECK IF COLUMN FIT (EXISTS) C CALL UTCFND(IDOUT,'FIT',1,COLIDS,ISTAT) IF(ISTAT.NE.0)THEN CALL UTCDEF(IDOUT,'FIT',' ',' ',TYREAL,1,COLIDS, * ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error adding column to '//TABLE GO TO 999 ENDIF ENDIF C C WRITE FITTED DATA C CALL UTCPTD(IDOUT,COLIDS(1),1,NS,FIT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to table '//TABLE GO TO 999 ENDIF CALL UTTCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing table '//TABLE GO TO 999 ENDIF ENDIF C C DONE C GO TO 1000 999 CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT) 1000 RETURN END