SUBROUTINE BSPEC C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C PROGRAM B_SPEC Fortran 77 MIDAS environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940331 0.0 M_SPEC C TH.MUELLER ST-ECF ESO GARCHING 940331 0.0 M_SPEC C M.R.ROSA ST-ECF ESO GARCHING 940504 1.0 B_SPEC C M.D. De La Pena STSCI 10 March 1998: changed DCOSD/DSIND to COS/SIN C with radian argument, use only COS/SIN. C M.D. De La Pena STSCI 01 March 2000: removed the use of IMPLICIT NONE. C C.PURPOSE C Model grating spectrograph for FOS scattered light analysis C Simplified code using blaze-function composition instead of C the full BLAZE*INTERFERENCE model of M_SPEC. Also simplified C LSF (Voigt profile + flat component) C Compute several items and combinations thereof C C Option 1: Intensity distribution input spectrum --> blazed grating C 2: Aperture diffraction pattern C 3: Convolution of blazed grating with aperture C 4: Blaze function C 0: option 1+3 C C.INPUTS C MIDAS character keywords IN_I and IN_L provide table names C MIDAS table "IN_I" provides input from columns C :angle = array of diffracted angles [deg] C :blg = array of intensity due to blazed grating only C (preferably [cts/s/angle step]) C :apd = array of intensity due to finite aperture (LSF) C MIDAS table "IN_L" provides input from columns C :wave = array of wavelengths [A] C :flx = array of intensities in inputs spectrum C (preferably [cts/s/A]) C MIDAS keyword RBUFF(14) provides parameters C 1> gsep = grating groove separatin in [mu] C 2> ntot = effective number of grooves C 3> blza = blaze angle [deg] C 4> inca = incident angle [deg] C C Best adjusted values as of 05-05-94 (M.R.) C C FOS Blue gsep ntot blza inca C G130H 1.000 34200 3.78 7.766 C G190H 1.458 23400 3.75 7.716 C G270H 2.083 16400 3.65 7.643 C G400H 3.055 11200 3.80 7.610 C G570H 4.444 7700 3.80 7.430 C C FOS Red gsep ntot blza inca C G190H 1.458 23400 3.74 7.684 C G270H 2.083 16400 3.65 7.649 C G400H 3.055 11200 3.80 7.607 C G570H 4.444 7700 3.80 7.426 C G780H 5.834 4900 3.80 7.701 C C 5> nlord= lowest order to be included (1.....) C 6> nuord= highest order to be included (...5...); 5 is typical C C 7> glsf = characteristic fall off of Voigt LSF (diode space) C 8> slsf = characteristic width of Voigt LSF (diode space) C 9> flat = flat component of LSF to simulate dust+roughness C C Typical values as of 05-05-94 C glsf slsf flat C 1" round apt: 0.003 0.01 1.5E-06 C C (if -9999. [deg] mapping will be centered) C 10> lmin = lower limit to input flux range of interest C default = :wave(1) C 11> lmax = upper limit to input flux range of interest C default = :wave(nrow(IN_L table)) C 12> iop = operation code C 0: evaluate grating, lsf and convolve C 1: evaluate grating scaled with flux --> :blg C 2: evaluate grating+aperture lsf --> :apd C 3: convolve :blg with :apd --> :ablg C 4: evaluate blaze funct of order nlord --> :bfu C C.OUTPUT C MIDAS table "IN_I" receives output into columns C :blg = array of intensity due to blazed grating only C (preferably [cts/s/angle step]) C :apd = array of intensity due to finite aperture C :ablg = array of intensity due to blazed grating + aperture C (preferably [cts/s/angle step]) C :benv = input spectrum dispersed with grating equation 1st ord C :bfu = blaze function (of selected order) C :gwav = wavelengths corresponding to diffracted angles 1st ord C C MIDAS table "IN_I" receives output into R*4 descriptor GRATING C copy of RBUFF to record parameters used C C-------------------------------------------------------------------------- C C ********************* C * Programmers Note: * C ********************* C C This MAIN is the interface to MIDAS keywords and MIDAS table environment C C 1) Access 2 tables, one containing at least 7 R*8 columns used as in/out, C the other one containing 2 R*4 columns, only input C C a) Typically the tables will have 2**7 - 2**16 rows C b) The number of elements (rows) need not necessarily be powers of 2 C because the convolution routine inserts into next fitting power2 C arrays. C C 2) Get 12 parameters into RBUFF C 3) Get VMR DOUBLE PREC storage for 2 phase arrays C 4) Hand over all that to the driving subroutine SPECTRUM C C SUBR. SPECTRUM and all SUBR. and FUNCTIONS called are C free of environment dependencies,with the following C --> EXCEPTION C SUBR. MCONVOL requires additional DOUBLE PREC. storage, C which is allocated using virtual memory mapping into C the COMMON /VMR/MEMD through SUBR. UDMGET C C 5) For other environments, eg. IRAF/STSDAS C a) Replace environmentals in MAIN to fit data structures C b) Replace the environment call in the VMR mapping routine C UDMGET by equivalents of STFXMP C C-------------------------------------------------------------------------- C INTEGER PH1PNT,PH2PNT INTEGER I,ISTAT,KNULL INTEGER NIR,NLR INTEGER ITID,ICI(7),ADI(7),LTID,ICL(2),ADL(2) INTEGER ITYP(7) REAL*4 RBUFF(12) REAL*4 DELTA CHARACTER*60 ITAB,LTAB CHARACTER*5 ICOL(7),LCOL(2),IFRM(7) CHARACTER*8 RLAB(12) CHARACTER*9 ICLU(7) CHARACTER*150 CONTXT CHARACTER*4 DET CHARACTER*5 FGWA CHARACTER*64 PST C C------------------------------------------------------------------------------ C Get IRAF MEM common into main program. C LOGICAL MEMB(1) INTEGER*2 MEMS(1) INTEGER*4 MEMI(1) INTEGER*4 MEML(1) REAL MEMR(1) DOUBLE PRECISION MEMD(1) COMPLEX MEMX(1) EQUIVALENCE (MEMB, MEMS, MEMI, MEML, MEMR, MEMD, MEMX) COMMON /MEM/ MEMD C C Codes for data types C INTEGER TYBOOL, TYREAL, TYDOUB PARAMETER (TYBOOL = 1) PARAMETER (TYREAL = 6) PARAMETER (TYDOUB = 7) C C Parameters relevant to table I/O C INTEGER TBALLR PARAMETER (TBALLR = 3) INTEGER TBMXCL PARAMETER (TBMXCL = 7) C C UMSPUT DESTINATIONS C INTEGER STDOUT PARAMETER (STDOUT = 1) INTEGER STDERR PARAMETER (STDERR = 2) C C------------------------------------------------------------------------------ C DATA ICOL/'angle','blg','apd','ablg','benv','bfu','gwav'/ DATA ICLU/'degrees',' ',' ',' ',' ',' ','angstroms'/ DATA IFRM/7*'E13.6'/ DATA ITYP/7*TYREAL/ DATA LCOL/'wave','flx'/ DATA RLAB/'gsep','ntot','blza','inca','low_ord','high_ord', + 'glsf','slsf','flat','lmin','lmax','iop'/ C C get table names and parameters from CL C call uclgst ('spectab',LTAB,ISTAT) call uclgst ('scattab',ITAB,ISTAT) call uclgst ('detector',det,ISTAT) call uclgst ('grating',fgwa,ISTAT) call uclgst ('instpars',pst,ISTAT) if(pst.eq.' ')pst='h'//fgwa(2:3)//det(1:1) call getpst(pst,rbuff(1),rbuff(2),rbuff(3),rbuff(4)) DO I = 5, 12 CALL UCLGSR (RLAB(I),RBUFF(I),ISTAT) END DO C C access the tables and map columns into virtual memory C CALL UTTINN (ITAB,ITID,ISTAT) NIR = 2048 CALL UTPPTI (ITID,TBALLR,NIR,ISTAT) CALL UTPPTI (ITID,TBMXCL,7,ISTAT) CALL UTCDEF (ITID,ICOL,ICLU,IFRM,ITYP,7,ICI,ISTAT) CALL UTTCRE (ITID,ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR creating table ' // ITAB GO TO 9999 ENDIF C DO I = 1, 7 CALL UDMGET (NIR, TYDOUB, ADI(I), ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR allocating memory for ADI' GO TO 9999 ENDIF ENDDO C CALL UTTOPN (LTAB,2,LTID,ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR opening table ' // LTAB GO TO 9999 ENDIF C CALL UTPGTI (LTID,21,NLR,ISTAT) CALL UDMGET (NLR, TYBOOL, KNULL, ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR allocating memory for KNULL' GO TO 9999 ENDIF C CALL UTCFND (LTID,'wave',1,ICL(1),ISTAT) IF (ICL(1).LE.0) CALL UTCFND (LTID,'wavelength',1,ICL(1),ISTAT) IF (ICL(1).LE.0) THEN CONTXT = 'ERROR locating wavelength column in table ' + // LTAB GOTO 9999 ELSE CALL UDMGET (NLR, TYREAL, ADL(1), ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR allocating memory for ADL' GO TO 9999 ENDIF CALL UTCGTR (LTID,ICL(1),1,NLR,MEMR(ADL(1)),MEMB(KNULL), + ISTAT) ENDIF C CALL UTCFND (LTID,'flx',1,ICL(2),ISTAT) IF (ICL(2).LE.0) CALL UTCFND (LTID,'flux',1,ICL(2),ISTAT) IF (ICL(2).LE.0) THEN CONTXT = 'ERROR locating flux column in table ' + // LTAB GOTO 9999 ELSE CALL UDMGET (NLR, TYREAL, ADL(2), ISTAT) IF( ISTAT .NE. 0) THEN CONTXT = 'ERROR allocating memory for ADL' GO TO 9999 ENDIF CALL UTCGTR (LTID,ICL(2),1,NLR,MEMR(ADL(2)),MEMB(KNULL), + ISTAT) ENDIF CALL UDMFRE (KNULL, TYBOOL, ISTAT) C C get dynamic storage for the phase arrays C CALL UDMGET(NIR,TYDOUB,PH1PNT,ISTAT) CALL UDMGET(NIR,TYDOUB,PH2PNT,ISTAT) C C set up diffracted angle array C DELTA = 45.0 / (NIR-1) DO I = 1, NIR MEMD(ADI(1)+I-1) = DBLE( I*DELTA - 15.0 ) END DO C C hand over to driving subroutine C CALL SPECT2(NIR,MEMD(ADI(1)),MEMD(ADI(2)), 1 MEMD(ADI(3)),MEMD(ADI(4)), 2 MEMD(ADI(5)),MEMD(ADI(6)), 3 MEMD(ADI(7)), 4 MEMD(PH1PNT),MEMD(PH2PNT), 5 NLR,MEMR(ADL(1)),MEMR(ADL(2)),RBUFF) C C write computed values to output columns C do i=1,7 call utcptd(itid,ici(i),1,nir,memd(adi(i)),istat) end do C C write descriptor to record parameters C DO I = 1, 12 CALL UTHADR (ITID,RLAB(I),RBUFF(I),ISTAT) END DO C GO TO 10000 C 9999 CALL UMSPUT( CONTXT, STDOUT+STDERR, 0, ISTAT) C C close the tables C 10000 CALL UTTCLO (ITID,ISTAT) CALL UTTCLO (LTID,ISTAT) C C free dynamic memory C DO I = 1, 7 IF (ADI(I).NE.0) CALL UDMFRE (ADI(I), TYDOUB, ISTAT) ENDDO DO I = 1, 2 IF (ADL(I).NE.0) CALL UDMFRE (ADL(I), TYREAL, ISTAT) ENDDO IF (PH1PNT.NE.0) CALL UDMFRE(PH1PNT,TYDOUB,ISTAT) IF (PH2PNT.NE.0) CALL UDMFRE(PH2PNT,TYDOUB,ISTAT) C END SUBROUTINE SPECT2(NDA,DA,BLGSPEC,APDA,ABLG, 1 BLGENV,BFU,GWAV,PH1,PH2,NW,WAVE,FLUX,RBUFF) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C SUBROUTINE SPECT2 Fortran 77 no environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940310 1.0 C TH.MUELLER ST-ECF ESO GARCHING 940310 1.0 C M.R.ROSA ST-ECF ESO GARCHING 940504 2.0 B_SPEC C C.PURPOSE C Model grating spectrograph for FOS scattered light analysis C Fast version of MODSPEC.FOR (DOUBLE PRECISION) C C Calls for diffrected angles in INPUT array DA(NDA) C subroutine BLAGRA to evaluate blazed grating intensities C subroutine LSFUNC to evaluate aperture diffraction function C subroutine MCONVOL to convolve intermediate results C subroutine BLAFU2 to evaluate blazed function C C.NOTES C USEAGE: C For studies of the FOS one might want to evaluate the two C components GRATING, APERTURE and their CONVOLUTION separately. C It is also useful to run the CONVOLUTION part in isolation on C the GRATING INTENSITY evaluated previously combined with an C LSF that has been determined by another program and fed into C column :apd of the input table. C C For production runs the program may be fed from a more elaborate C command procedure or main module that have a detailed knowledge C of FOS components and the yet to be made experiences about LSFs. C C PREPARATION of FLUX data: C The program is intended to produce output data that correspond to C the (dead-time and dead-diode corrected) COUNT spectra (.c4h) C of FOS data sets. The working independent is diffracted angle [deg], C which is a linear function of the diode [1-512] resp. detector pixel C space (1 diode = 0.0058 [deg], diode(1) = -1.47[deg]. Diffracted C angle 0[deg] corresponds to diode [256]. C C This requires that the input FLUX is already scaled by all optical C components other than the gratings blaze function, in particular C the conversion from F(lambda)d(lambda) to cts/pix has to made C before entering. Assuming the input spectrum is F(l)dl this has to C be multiplied by the reflectance of OTA, GRAZING MIRROR, COLLIMATOR, C the DETECTOR DQE, an undispersed light reflectance characteristic C for the grating (the COLLIMATOR seems to be a good template) and C an aperture throughput function (light loss at colimator due to C diffraction and light loss at aperture due to cutting the HST PSF). C NOTE THAT ALL THIS COMBINED IS NOT THE IVS CURVE, since the ivs C curve in addition compensates for the blaze function. ALSO - one C needs to get all wavelengths onto the grating, properly weighted, - C while the ivs truncates the light for wavelengths outside the C range seen by the detector. The blaze function can be studied using C option 4 to call BLAFU2. C C SAMPLING: C The independent array DA (diffracted angle) should have C sufficient sampling to oversample the resolution at the detector. C For FOS 1 diode corresponds to 0.0058 [deg] and the detector C spans -1.47 to +1.47 [deg]. Zero's order is at -incangle = about C 7.5-7.8 [deg] (high res), or -0.5 [deg] for low res gratings. C It is adviseable to cover at least the range -9 to +15 [deg] in C order to get the contribution of orders 0,1,2 to scattered light. C C The input flux data need not be highly oversampled. The LSFs core C has a FWHM of at least 1 diode = 0.0058 [deg] even for the smallest C apertures. To evaluate the level of scattered light even a sampling C of 10-50 diodes is sufficient. C C C.INPUTS/OUTPUTS C NDA No of diffracted angles C DA diffracted angles[degree] C BLGSPEC diffracted intensity C PH1,PH2 phasedifferences C APDA grating+aperture lsf function C ABLG convolution of lsf and blgspec C BLGENV input flux dispersed with grating equation for first order C GWAV wavelengths corresponding to diffracted angles in first order C BFU blaze function C NW No of wavelength lambda C WAVE wavelengths lambda[A] C FLUX flux at lambda C RBUFF grating parameters and indices and options C 1=s[mu] 2=ngrooves 3=blazang 4 incang C 5,6=order limits 7,8,9=voigt parameters g,s, and flat lsf component C angle 10,11=lower,upper lambda 12=program option C ----------------------------------------------------------------------- C INTEGER*4 I,J,NW,NDA,INDA,N1,N2,NDLO,NDUP,IOP,NTLO,NTUP INTEGER*4 NACT,NLORD,NUORD,INDZERO REAL*4 RBUFF(12),LAMB,WAVE(NW),FLUX(NW),ALINT2,FMIN, * FMAX DOUBLE PRECISION BLGSPEC(NDA),DA(NDA),PH1(NDA),PH2(NDA) DOUBLE PRECISION FLAM,PI,WMIN,WMAX,GSEP,NGROO,GS,TOTF DOUBLE PRECISION BLZANGLE,INCANGLE,SININCA,COSBLZA,SINTHETA DOUBLE PRECISION APDA(NDA),ABLG(NDA),CA,GLSF,SLSF,FLAT DOUBLE PRECISION BLGENV(NDA),BFU(NDA),GWAV(NDA),ORD,ZERO,TORD DOUBLE PRECISION DGTORD character*160 contxt c integer istat C EXTERNAL ALINT2 DATA PI /3.14159265358979323846D+00/ C DGTORD = PI / 180.0D0 C C Prepare parameters ! typical GSEP = DBLE(RBUFF(1))*1.0D-06 ! s in meters 2.E-6 GS = GSEP*1.D+10 ! s in Angstrom 20000 NGROO = DBLE(RBUFF(2)) ! # of grooves 10000 BLZANGLE = DBLE(RBUFF(3)) ! blaze angle 4 deg INCANGLE = DBLE(RBUFF(4)) ! incident angle 7.5 deg NLORD = INT(RBUFF(5)) ! low order 1 NUORD = INT(RBUFF(6)) ! high order 5 WMIN = DBLE(MAX(WAVE(1),RBUFF(10))) ! wavelength - 1000 A WMAX = DBLE(MIN(WAVE(NW),RBUFF(11))) ! limts for flux 9000 A IOP = INT(RBUFF(12)) ! Options 0..4 c TYPE *,RBUFF write(contxt,10) rbuff 10 format(12e13.5) c call umsput(contxt,1,0,istat) SINTHETA = SIN(DGTORD*(INCANGLE-BLZANGLE)) SININCA = SIN(DGTORD*INCANGLE) COSBLZA = COS(DGTORD*BLZANGLE) C C Fill phase and lambda arrays DO I = 1,NDA PH1(I) = COSBLZA*(SIN(DGTORD*(DA(I)-BLZANGLE))+SINTHETA) ! slow phase PH2(I) = SIN(DGTORD*DA(I))+SININCA ! fast phase GWAV(I) = GS*PH2(I) ! wavelengths 1st ord ENDDO C C IOP = 0 produces 1,2 and 3 IF(IOP.EQ.1) GOTO 1000 ! --> grating only IF(IOP.EQ.2) GOTO 2000 ! --> aperture only IF(IOP.EQ.3) GOTO 3000 ! --> convolution only IF(IOP.EQ.4) GOTO 4000 ! --> blaze function only C C Prepare output arrays and indices for grating evaluation 1000 NDLO = 1 NDUP = NDA NTLO = NDA NTUP = 1 DO I = 1,NDA BLGSPEC(I) = 0.0D0 BLGENV(I) = 0.0D0 IF ((GS*PH2(I)).LT.WMIN) NDLO = I ! index for 1st order wavelength min IF ((GS*PH2(I)).LE.WMAX) NDUP = I ! index for 1st order wavelength max IF ((GS*PH2(I)).LE.(-WMIN)) NTLO = I ! index for -1 ord wavelength min IF ((GS*PH2(I)).LT.(-WMAX)) NTUP = I ! index for -1 ord wavelength max IF (DA(I).LT.0.0D0) INDA=I ! index for da = 0. IF (DA(I).LE.(-INCANGLE)) INDZERO = I ! index for zero order ENDDO c TYPE *,INDZERO,DA(INDZERO) write(contxt,15)indzero,da(indzero) 15 format(i5,f8.3) c call umsput(contxt,1,0,istat) c TYPE *,'angle indices: ',NDLO,NDUP,' wavelength: ',WMIN,WMAX write(contxt,20)ndlo,ndup,wmin,wmax 20 format('angle indices: ',2i5,' wavelength: ',2f10.3) call umsput(contxt,1,0,istat) C C Fill the 1st and -1st order dispersed input flux (blaze not included) C NACT = 1 !start index for flux interpolation N2 = NW !end index DO I = NDLO,NDUP N1 = NACT LAMB = SNGL(GS*PH2(I)) ![A] BLGENV(I) = DBLE(ALINT2(NW,WAVE,FLUX,LAMB,N1,N2,NACT)) ENDDO NACT = 1 !start index for flux interpolation N2 = NW !end index DO I = NTLO,NTUP,-1 N1 = NACT LAMB = -SNGL(GS*PH2(I)) ![A] BLGENV(I) = DBLE(ALINT2(NW,WAVE,FLUX,LAMB,N1,N2,NACT)) ENDDO CALL TOTAL(BLGENV,NDA,1,TOTF,FMIN,FMAX) C C Go through the Flux(angle(lambda)) space, evaluate the blaze function C function at each diffracted angle, multiply with the interpolated C Flux(lambda) for each order (note dilution by lambda-spread) and sum up. C Lambda(order) = gsep*(sin(diffangle)+sin(incangle))/order C Take care of decreasing lambda at diffracted angles shortward of 0 order C CALL BLAFU2(NDA,DA,DBLE(RBUFF(4)),DBLE(RBUFF(3)),BFU,1.0D0) C DO J = NLORD,NUORD ORD = 1.0D0*J NACT = 1 N2 = NW DO I = INDZERO,NDA N1 = NACT LAMB = SNGL(GS*PH2(I)/ORD) ![A] FLAM = DBLE(ALINT2(NW,WAVE,FLUX,LAMB,N1,N2,NACT)) BLGSPEC(I) = BLGSPEC(I)+FLAM*BFU(I)/ORD ENDDO NACT = 1 N2 = NW DO I = INDZERO-1,1,-1 N1 = NACT LAMB = -SNGL(GS*PH2(I)/ORD) ![A] FLAM = DBLE(ALINT2(NW,WAVE,FLUX,LAMB,N1,N2,NACT)) BLGSPEC(I) = BLGSPEC(I)+FLAM*BFU(I)/ORD ENDDO c TYPE *,'Computing order ',J write(contxt,25)j 25 format('Computing order ',i3) call umsput(contxt,1,0,istat) ENDDO C C Zero order flux is (input - flux dispersed into non-zero orders). C Usually considering orders 1 to 5 (and -1 to -5) is sufficient C ZERO = 0.0D0 TORD = 0.0D0 DO I = 1,NDA ZERO = ZERO+(BLGENV(I)-BLGSPEC(I)) TORD = TORD+BLGSPEC(I) ENDDO BLGSPEC(INDZERO) = ZERO c TYPE *,'Flux in zero order:',BLGSPEC(INDZERO),' in higher orders:',TORD write(contxt,30)blgspec(indzero),tord 30 format('Flux in zero order:',e13.5,' in higher orders:',e13.5) call umsput(contxt,1,0,istat) C C C This was the grating part C Now evaluate the aperture diffraction LSF centered in DA array, C shift to angle CA if required for inspection. You better keep C it centered for proper results during convolution. C 2000 IF(IOP.GT.0.and.IOP.NE.2) GOTO 3000 GLSF = DBLE(RBUFF(7)) SLSF = DBLE(RBUFF(8)) FLAT = DBLE(RBUFF(9)) CA = DA(INT(NDA/2+0.5)) CALL LSFUNC(NDA,DA,APDA,CA,GLSF,SLSF,FLAT) C C Convolve the grating part with the aperture LSF. C If LSF was not correctly centered -- too bad C 3000 IF(IOP.GT.0.AND.IOP.NE.3) GOTO 4000 CALL MCONVL(NDA,BLGSPEC,APDA,ABLG) C C Evaluate the blaze function = intensity distribution(lambda) C 4000 IF(IOP.NE.4) RETURN ORD = 1.0D0*NLORD CALL BLAFU2(NDA,DA,DBLE(RBUFF(4)),DBLE(RBUFF(3)),BFU,ORD) C RETURN END C REAL FUNCTION ALINT2(N,X,Y,XC,N1,N2,IAL) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C SUBROUTINE ALINT2 Fortran 77 no environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940310 1.0 C TH.MUELLER ST-ECF ESO GARCHING 940310 1.0 C C.PURPOSE C Linear interpolation in logarithmic data space in REAL*4 C Assumes continously increasing indepent, handles search limits, C returns actual lower index for calls from loops C C.INPUTS C N dimension of indepent,dependent data array C X,Y indepent,dependent data arrays C XC value of independent at which to evaluate C N1,N2 indices of search limits C C.OUTPUTS C ALINT2 interpolated value C IAL actual lower index limit in independent C ------------------------------------------------------------------- INTEGER*4 N,N1,N2,IAL,I REAL*4 X(N),Y(N),XC,Y1,Y2 DO 1 I=N1,N2 1 IF(X(I).GT.XC) GOTO 2 2 ALINT2 = 1.E-29 IF(Y(I-1).LE.0.0.OR.Y(I).LE.0.0) RETURN Y1=ALOG10(Y(I-1)) Y2=ALOG10(Y(I)) ALINT2=10.**(Y1+(Y2-Y1)*(XC-X(I-1))/(X(I)-X(I-1))) IAL = MAX(I-1,1) RETURN END C SUBROUTINE BLAFU2(NDA,DA,INCA,BLZA,BFU,ORD) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.MODULE C SUBROUTINE BLAFU2 Fortran 77 no environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940310 1.0 C TH.MUELLER ST-ECF ESO GARCHING 940310 1.0 C M.R.ROSA ST-ECF ESO GARCHING 940504 2.0 B_SPEC C C.PURPOSE C Model grating for FOS scattered light analysis C C Evaluate the blaze function of a reflection grating C at wavelengths corresponding to diffracted angles DA C C bfu = [sinc(g)]**2 C g = pi * cos(d) * (sin(b-d)+sin(a-d)) / (sin(b)+sin(a)) C a = incident angle C b = diffracted angle C d = blaze angle C C.INPUTS/OUTPUTS C NDA No of diffracted angles C DA diffracted angles [degree] C INCA incident angle C BLZA blaze angle C BFU blaze function C ORD order C ------------------------------------------------------------------- INTEGER*4 NDA,I DOUBLE PRECISION DA(NDA),BFU(NDA),ORD DOUBLE PRECISION PI,GAMMA,DIVI,PBLZ,SINIB,SINI DOUBLE PRECISION BLZA,INCA,INTS,BFUNC DOUBLE PRECISION DGTORD DATA PI /3.14159265358979323846D+00/ C DGTORD = PI / 180.0D0 C PBLZ = DABS(ORD)*PI*COS(DGTORD*BLZA) SINIB = SIN(DGTORD*(INCA-BLZA)) SINI = SIN(DGTORD*INCA) INTS = 1.0D0 DO I = 1,NDA DIVI = SIN(DGTORD*DA(I))+SINI IF (DABS(DIVI).LT.1.0D-38) DIVI = 1.0D0 GAMMA = PBLZ*(SIN(DGTORD*(DA(I)-BLZA))+SINIB)/DIVI CALL SINCF2(GAMMA,GAMMA,INTS,BFUNC) BFU(I) = BFUNC END DO C RETURN END C SUBROUTINE LSFUNC(NDA,DA,LSF,CA,GX,SX,FX) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C SUBROUTINE LSFUNC Fortran 77 no environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940504 1.0 B_SPEC C C.PURPOSE C Model aperture for FOS scattered light analysis C Voigt profile + flat background C C.INPUTS/OUTPUTS C NDA No of diffracted angles C DA diffracted angles[degree] C APDA aperture function C CA central angle C LAM wavelength C APSIZ aperture size C C ------------------------------------------------------------------- C C T.Mueller 28-Feb-94 C M.Rosa 15-Mar-94 C----------------------------------------------------------------------- INTEGER*4 NDA,I REAL*4 FMIN,FMAX DOUBLE PRECISION DA(NDA),LSF(NDA),CA,GX,SX,FX DOUBLE PRECISION PI,LN2,NORM,VLS,GD,SD,DX,DD,SDD C DATA PI /3.14159265358979323846D+00/ C c Type *,' LSF' c call umsput(' LSF',1,0,istat) C LN2 = DLOG(2.0D0) GD = GX*0.580D-02 ! * degrees per diode SD = SX*0.580D-02 ! * degrees per diode DX = GX/2.0D0+DSQRT(GX*GX/4.0D0+2.0D0*LN2*SX*SX) C DO I = 1,NDA DD = (DA(I)-CA)/DX SDD = DD*DD VLS = (1.0D0-GX/DX)*DEXP(-LN2*SDD)+GX/DX/(1.0D0+SDD) LSF(I) = VLS+FX ENDDO CALL TOTAL(LSF,NDA,1,NORM,FMIN,FMAX) CALL NORMAL(LSF,NDA,1,NORM) C RETURN END C SUBROUTINE SINCF2(DA,DB,ILA,SCF) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C SUBROUTINE SINCF2 Fortran 77 no environment C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940310 1.0 C TH.MUELLER ST-ECF ESO GARCHING 940310 1.0 C C.PURPOSE C compute DBLE precision SINC function C C.INPUTS/OUTPUTS C DA rad 1 C DB rad 2 C SINCF2 = [sin(da)/db]**2 C -------------------------------------------------------- DOUBLE PRECISION DA,DB,ILA,SCF IF (DABS(DB).LT.1.0D-10) THEN SCF = 1.0D0*ILA ELSE SCF = ((SIN(DA)/DB)**2)*ILA END IF RETURN END C SUBROUTINE MCONVL (NP,INPUT,PSF,RESULT) C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C SUBROUTINE MCONVOL Fortran 77 MIDAS VMR mapping req. C C.AUTHORS C M.R.ROSA ST-ECF ESO GARCHING 940310 1.0 C TH.MUELLER ST-ECF ESO GARCHING 940310 1.0 C R.HOOK ST-ECF ESO GARCHING 92-93 0.0 C C.PURPOSE C Convolve 2 DBLE data arrays (1D), dimension not necessarily C power of 2 C Uses R. Hook's code for Lucy deconv. in 2D C C Environmental dependency: COMMON /VMR/MEMD C CALL UDMGET C C.INPUTS/OUTPUTS C NP # of data points in the input/output arrays C INPUT DBLE input data array C PSF DBLE input psf array C RESULT DBLE output array of convolved data C C ------------------------------------------------------------------- c INTEGER MEMD(1) c COMMON /VMR/MEMD REAL*4 FMIN,FMAX INTEGER ISTAT,NX,NY,DIMS(2),I,NP,NPOW,NTEST INTEGER DATPNT,PSFPNT,RESPNT,WORK,PSFFFT DOUBLE PRECISION INPUT(NP),PSF(NP),RESULT(NP),PT,DT character*80 contxt C------------------------------------------------------------------------------ C Get IRAF MEM common into main program. C LOGICAL MEMB(1) INTEGER*2 MEMS(1) INTEGER*4 MEMI(1) INTEGER*4 MEML(1) REAL MEMR(1) DOUBLE PRECISION MEMD(1) COMPLEX MEMX(1) EQUIVALENCE (MEMB, MEMS, MEMI, MEML, MEMR, MEMD, MEMX) COMMON /MEM/ MEMD C------------------------------------------------------------------------------ C C First determine power2 array that fits data DO I=1,20 NTEST = INT(2.**I) IF(NTEST.LT.NP) NPOW = I ENDDO NX=2**(NPOW+1) NY=1 DIMS(1) = NX DIMS(2) = NY C Sorry - who dares to handle arrays larger 2**20 = 1.05E+06 C C Insert DATA and PSF arrays into power2 arrays for FFTing CALL UDMGET(NX*NY,7,DATPNT,ISTAT) CALL UINSAR(NP,INPUT,NX,NY,MEMD(DATPNT),ISTAT) CALL UDMGET(NX*NY,7,PSFPNT,ISTAT) CALL UINSAR(NP,PSF,NX,NY,MEMD(PSFPNT),ISTAT) C C Prepare the arrays for the Numerical Recipes FFT routine C and do the FFTs of the PSF and rotated PSF C C At this point we also need to check that the PSFs are normalised to C a total of 1 - if they aren't we divide C To avoid numerical trouble we do this also for the DATA, but will later C renormalize the RESULT CALL TOTAL(MEMD(DATPNT),NX,NY,DT,FMIN,FMAX) IF(DABS(DT-1.0D0).GT.1.0E-5) THEN c Type *,'DATA will be normalized:', DT write(contxt,10) dt 10 format('DATA will be normalized: ', d11.4) call umsput(contxt,1,0,istat) CALL NORMAL(MEMD(DATPNT),NX,NY,DT) ELSE c Type *,'DATA normal:', DT write(contxt,20) dt 20 format('DATA normal: ', d11.4) call umsput(contxt,1,0,istat) ENDIF C CALL TOTAL(MEMD(PSFPNT),NX,NY,PT,FMIN,FMAX) IF(DABS(PT-1.0D0).GT.1.0E-5) THEN c Type *,'PSF will be normalized:', PT write(contxt,30) pt 30 format('PSF will be normalized: ', d11.4) call umsput(contxt,1,0,istat) CALL NORMAL(MEMD(PSFPNT),NX,NY,PT) ELSE c Type *,'PSF normal:', PT write(contxt,40) pt 40 format('PSF normal: ', d11.4) call umsput(contxt,1,0,istat) ENDIF C C Now we need additional memory while FFTing CALL UDMGET(NX*NY*2,7,PSFFFT,ISTAT) CALL UDMGET(NX*NY*2,7,WORK,ISTAT) CALL UDMGET(NX*NY,7,RESPNT,ISTAT) C C We need to fill the PSFs FFT space CALL DFILL(MEMD(PSFPNT),NX,NY,MEMD(PSFFFT)) CALL DFOURN(MEMD(PSFFFT),DIMS,2,+1) C C Convolve the data with the PSF CALL DCONV(MEMD(DATPNT),NX,NY,MEMD(WORK), 1 MEMD(PSFFFT),MEMD(RESPNT),1) C C Extract the result and renormalize CALL UEXTAR(NX,NY,MEMD(RESPNT),NP,RESULT,ISTAT) DO I = 1,NP RESULT(I) = RESULT(I)*DT ENDDO C Done RETURN END C SUBROUTINE NORMAL(DATA,NX,NY,TOT) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Normalize DBLE array C----------------------------------------------------------------------------- INTEGER NX,NY,I,J DOUBLE PRECISION DATA(NX,NY),TOT C DO J=1,NY DO I=1,NX DATA(I,J) = DATA(I,J)/TOT ENDDO ENDDO C RETURN END C SUBROUTINE TOTAL(DATA,NX,NY,TOT,FMIN,FMAX) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Add up all elements in DBLE array, return total and limits C----------------------------------------------------------------------------- INTEGER NX,NY,I,J REAL*4 FMIN,FMAX DOUBLE PRECISION DATA(NX,NY),TOT,DMIN,DMAX C TOT=0.0D0 DMIN = 0.0D0 DMAX = 0.0D0 DO J=1,NY DO I=1,NX IF(DATA(I,J).LT.DMIN) THEN DMIN = DATA(I,J) ELSE IF(DATA(I,J).GT.DMAX) DMAX = DATA(I,J) ENDIF TOT=TOT+DATA(I,J) ENDDO ENDDO C FMIN = DMIN FMAX = DMAX C RETURN END C SUBROUTINE UINSAR(ND,DATA,NX,NY,ARR,ISTAT) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Insert DBLE data centered into an array for use in FFT convolution. C Make sure first and last pix is 0.0, pix 2 interpolated to pix 3, C pix NX-1 interpolated to pix NX-2. C---------------------------------------------------------------------------- INTEGER ISTAT,NX,NY,ND,NS,NE,I,J DOUBLE PRECISION ARR(NX,NY),DATA(ND) ISTAT = 0 C NS = INT(NX/2)-INT(ND/2)+1 NE = NS+ND C DO I = 1,MAX(NS-1,1) ARR(I,1) = 0.0D0 ENDDO DO I = MIN(NE+1,NX),NX ARR(I,1) = 0.0D0 ENDDO J = 0 DO I = NS,MIN(NE,NX) J = J+1 ARR(I,1) = DATA(J) ENDDO ARR(1,1) = 0.0D0 ARR(NX,1) = 0.0D0 ARR(2,1) = 0.5*ARR(3,1) ARR(NX-1,1) = 0.5*ARR(NX-2,1) C RETURN END C SUBROUTINE DFILL(DATA,NX,NY,OUT) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Fill up an array in the way required by the Numerical Recipes C FFT routines C C Note this routine has a 'fudge factor' of 1 pixel to align coordinate C systems correctly. C---------------------------------------------------------------------------- C INTEGER NX,NY,IP,I,J DOUBLE PRECISION DATA(NX,NY) DOUBLE PRECISION OUT(NX*NY*2) C IP=1 IF(NY.GT.1) THEN DO I=1,NX OUT(IP)=0.0D0 IP=IP+2 ENDDO DO J=1,NY-1 OUT(IP)=0.0D0 IP=IP+2 DO I=1,NX-1 OUT(IP-1)=0.0D0 OUT(IP)=DATA(I,J) IP=IP+2 ENDDO ENDDO ELSE DO I=1,NX OUT(IP)=DATA(I,1) OUT(IP+1)=0.0D0 IP=IP+2 ENDDO ENDIF C RETURN END C SUBROUTINE DFOURN(DATA,NN,NDIM,ISIGN) C C Numerical Recipes FFT routine - modified to use entirely C DOUBLE PRECISION throughout. C INTEGER NN,NDIM,ISIGN,NTOT,IDIM,NPREV,N,IP3,IP2,NREM INTEGER I2REV,I2,IP1,I1,I3,I3REV,IBIT,IFP1,IFP2,K1,K2 REAL*8 WR,WI,WPR,WPI,WTEMP,THETA,DATA,TEMPR,TEMPI DIMENSION NN(NDIM),DATA(*) C NTOT=1 DO 11 IDIM=1,NDIM NTOT=NTOT*NN(IDIM) 11 CONTINUE NPREV=1 DO 18 IDIM=1,NDIM N=NN(IDIM) NREM=NTOT/(N*NPREV) IP1=2*NPREV IP2=IP1*N IP3=IP2*NREM I2REV=1 DO 14 I2=1,IP2,IP1 IF(I2.LT.I2REV)THEN DO 13 I1=I2,I2+IP1-2,2 DO 12 I3=I1,IP3,IP2 I3REV=I2REV+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(I3REV) DATA(I3+1)=DATA(I3REV+1) DATA(I3REV)=TEMPR DATA(I3REV+1)=TEMPI 12 CONTINUE 13 CONTINUE ENDIF IBIT=IP2/2 1 IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN I2REV=I2REV-IBIT IBIT=IBIT/2 GO TO 1 ENDIF I2REV=I2REV+IBIT 14 CONTINUE IFP1=IP1 2 IF(IFP1.LT.IP2)THEN IFP2=2*IFP1 THETA=ISIGN*6.28318530717959D0/(IFP2/IP1) WPR=-2.D0*SIN(0.5D0*THETA)**2 WPI=SIN(THETA) WR=1.D0 WI=0.D0 DO 17 I3=1,IFP1,IP1 DO 16 I1=I3,I3+IP1-2,2 DO 15 I2=I1,IP3,IFP2 K1=I2 K2=K1+IFP1 TEMPR=WR*DATA(K2)-WI*DATA(K2+1) TEMPI=WR*DATA(K2+1)+WI*DATA(K2) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR DATA(K1+1)=DATA(K1+1)+TEMPI 15 CONTINUE 16 CONTINUE WTEMP=WR WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTEMP*WPI+WI 17 CONTINUE IFP1=IFP2 GO TO 2 ENDIF NPREV=N*NPREV 18 CONTINUE C RETURN END C SUBROUTINE DCONV(DATA,NX,NY,WORK,PSFFFT,OUT,FLAG) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Convolve a one or two dimensional DBLE array with a PSF which is C supplied as a double precision FFT. Uses Numerical Recipes C C Richard Hook, ST-ECF, October 1991 C Double precision version, March 1992 C 1D version, August 1992 C---------------------------------------------------------------------------- C INTEGER I,J,N,IP,IQ DOUBLE PRECISION A,B,C,D,FAC INTEGER NX,NY ! ARRAY SIZE INTEGER NDIMS,DIMS(2),NX2,NY2 INTEGER FLAG ! ROTATED OR NOT? REAL*8 DATA(NX,NY) ! INPUT DATA ARRAY REAL*8 OUT(NX,NY) REAL*8 PSFFFT(NX*NY*2) ! FFT OF PSF REAL*8 WORK(NX*NY*2) ! WORKSPACE FOR FFT OF DATA REAL*8 DFLAG C DFLAG=DBLE(FLAG) IF(NY.EQ.1) THEN NDIMS=1 ELSE NDIMS=2 ENDIF C C Copy the data into the real part of the working array IP=1 DO J=1,NY DO I=1,NX WORK(IP)=DATA(I,J) WORK(IP+1)=0.0D0 IP=IP+2 ENDDO ENDDO C C Take the forward FFT of the data array DIMS(1)=NX DIMS(2)=NY C CALL DFOURN(WORK,DIMS,NDIMS,+1) C C Do a complex multiplication of the FFT of the data and of the PSF C C If the rotation flag is set (to -1) we multiply the imaginary part C of the FFT by -1 IP=1 IQ=2 DO N=1,NX*NY A=WORK(IP) B=WORK(IQ) C=PSFFFT(IP) D=DFLAG*PSFFFT(IQ) WORK(IP)=A*C-B*D WORK(IQ)=A*D+B*C IP=IP+2 IQ=IQ+2 ENDDO C C Take the inverse FFT CALL DFOURN(WORK,DIMS,NDIMS,-1) C C Copy the result into the output array C and 'quadrant swap': FAC=DFLOAT(NX*NY) IP=1 C NX2=MAX(1,NX/2) NY2=MAX(1,NY/2) C IF(NY.GT.1) THEN DO J=1,NY2 DO I=1,NX2 OUT(I+NX/2,J+NY/2)=WORK(IP)/FAC IP=IP+2 ENDDO IP=IP+NX ENDDO C IP=NX+1 DO J=1,NY2 DO I=1,NX2 OUT(I,J+NY/2)=WORK(IP)/FAC IP=IP+2 ENDDO IP=IP+NX ENDDO C IP=NX*NY+1 DO J=1,NY2 DO I=1,NX2 OUT(I+NX/2,J)=WORK(IP)/FAC IP=IP+2 ENDDO IP=IP+NX ENDDO C IP=NX*NY+NX+1 DO J=1,NY2 DO I=1,NX2 OUT(I,J)=WORK(IP)/FAC IP=IP+2 ENDDO IP=IP+NX ENDDO ELSE C C 1D case DO I=1,NX2 OUT(I+NX2,1)=WORK(IP)/FAC IP=IP+2 ENDDO C DO I=1,NX2 OUT(I,1)=WORK(IP)/FAC IP=IP+2 ENDDO ENDIF C RETURN END C SUBROUTINE UEXTAR(NX,NY,ARR,ND,DATA,ISTAT) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Extract data centered from an DP array C----------------------------------------------------------------------------- INTEGER ISTAT,NX,NY,ND,NS,NE,I,J DOUBLE PRECISION ARR(NX,NY),DATA(ND) ISTAT = 0 C NS = INT(NX/2)-INT(ND/2)+1 NE = NS+ND C J = 0 DO I = MAX(1,NS),MIN(NE,NX) J = J+1 DATA(J) = ARR(I,1) ENDDO C RETURN END