C @(#)spec3g.for 17.1.1.1 (ES0-DMD) 01/25/02 17:14:57 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 SPEC3G C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C C Program SPEC3G version 1.00 860908 C version 1.01 870609 C (Allow any no. of i/p lines.) C version 2.00 890424 (pMIDAS) C F. Murtagh STECF C M. Peron IPG 890901 add include files C P. Ballester IPG 910313 Compil. option -u C and more C C.KEYWORDS C Simulated images, test images. C C.PURPOSE C Create spectrum of specified lines, each a gaussian in form. C C.OUTPUT C C Keys: OUT_A/C/1/60 output data array C INPUTR/R/1/1 centring relative to pixel centre C IN_TAB/C/1/60 table of positions of lines C INPUTR/R/2/1 full width half max. C C----------------------------------------------------------- C C IMPLICIT NONE C REAL RMIN,RMAX REAL STEPO(3),STARTO(3),CUTS(4) DOUBLE PRECISION DSTEP(3),DSTART(3) INTEGER NPIXO(3),KUN,KNUL INTEGER MADRID CHARACTER*60 OUTIMA,INTAB CHARACTER*72 IDENTO CHARACTER*80 CUNITO C INTEGER IACT,ISTAT,NDIMO,IDIM INTEGER*8 JPNTR INTEGER IMNO,TID INTEGER IACTV,NC,NR,NS,NAC,NAR INTEGER*8 IPTR INTEGER KUNIT REAL CENTR,FWHM C INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C C C ... get into MIDAS C CALL STSPRO('SPEC3G') C C ... get name of output frame C CALL STKRDC('OUT_A',1,1,60,IACT,OUTIMA,KUN,KNUL,ISTAT) C C ... get centring. C CALL STKRDR('INPUTR',1,1,IACT,CENTR,KUN,KNUL,ISTAT) C C ... get gaussian fwhm. C CALL STKRDR('INPUTR',2,2,IACT,FWHM,KUN,KNUL,ISTAT) IF (FWHM.LT.0.2) THEN CALL STTPUT(' The half maximum is too small (< 0.2).',ISTAT) CALL STTPUT(' Aborting.',ISTAT) GOTO 1000 ENDIF C C ... map output image C NDIMO = 1 IDIM = 660 NPIXO(1) = IDIM STARTO(1) = 1.0 STEPO(1) = 1.0 CUNITO = ' NONE' IDENTO = ' ARTIFICIAL SPECTRUM' DSTART(1) = STARTO(1) DSTEP(1) = STEPO(1) CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NDIMO,NPIXO,DSTART,DSTEP, . IDENTO,CUNITO,JPNTR,IMNO,ISTAT) C C ... get line positions from table. C CALL STKRDC('IN_TAB',1,1,60,IACTV,INTAB,KUN,KNUL,ISTAT) CALL TBTOPN(INTAB,0,TID,ISTAT) CALL TBIGET(TID,NC,NR,NS,NAC,NAR,ISTAT) IF (NC.LT.1.OR.NR.LT.1) THEN CALL STTPUT(' Invalid nos. of rows or columns.',ISTAT) CALL STTPUT(' in table specifying lines.',ISTAT) CALL STTPUT(' Aborting.',ISTAT) GOTO 1000 ENDIF CALL TBCMAP(TID,1,IPTR,ISTAT) C C ... now do the work C CALL PATTERN(MADRID(JPNTR),RMIN,RMAX,NR,MADRID(IPTR),CENTR,FWHM, X ISTAT) C C ... write cuts C CUTS(1) = RMIN CUTS(2) = RMAX CUTS(3) = RMIN CUTS(4) = RMAX CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUNIT,ISTAT) C C ... end C CALL STFCLO(IMNO,ISTAT) 1000 CONTINUE CALL STSEPI END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C C Program PATTERN C F. MURTAGH ST-ECF Version 1.0 860411 C C.KEYWORDS C C Test patterns, simulated images. C C C.OUTPUT PARAMETERS C C ARR = the frame, C RMIN, RMAX = cut values (max. and min. flux values). C C---------------------------------------------------------------------- SUBROUTINE PATTERN(ARR,RMIN,RMAX,N,VALS,CENTR,FWHM, X ISTAT) C IMPLICIT NONE C REAL ARR(660),VALS(1) REAL GAUSS1(220),GAUSS2(220),GAUSS3(220) C INTEGER N,ISTAT,I,II,IVAL,III REAL RMIN,RMAX,CENTR,FWHM,ZL,ZU,Y,SD,VAL,SCAL C C C DO I = 1, 220 GAUSS1(I) = 0.0 ! Initialize GAUSS2(I) = 0.0 ! Initialize GAUSS3(I) = 0.0 ! Initialize ENDDO SD = FWHM/2.354 ! Det. std. dev. from half width. DO I = 1, 660 ARR(I) = 0.0 ! Initialize ENDDO C C ... Det. flux for a Gaussian, with mean at (of offset from) 110.0. C DO I = 1,220 ZL = FLOAT(I)-0.5 ! Pixel extremity (lower). ZU = FLOAT(I)+0.5 ! Pixel extremity (upper). Y = 110.0 ! Mean. CALL AG(ZL,ZU,Y,SD,VAL) ! Det. flux in given pixel. GAUSS1(I) = VAL Y = 110.5 ! Mean offset by 0.5 pix. unit CALL AG(ZL,ZU,Y,SD,VAL) GAUSS2(I) = VAL Y = 110.0+CENTR ! Mean offset by user-def'd. val. CALL AG(ZL,ZU,Y,SD,VAL) GAUSS3(I) = VAL ENDDO C C ... Do for segments 1-220,221-440,441-660 of ARR. C DO II = 1,N IVAL = VALS(II) C DO I = 1, 220 III = I-110+IVAL IF (III.LE.0) GOTO 100 IF (III.GT.220) GOTO 100 ARR(III) = ARR(III) + GAUSS1(I) ! Segment 1-200 ARR(III+220) = ARR(III+220) + GAUSS2(I) ! Segment 201-440 ARR(III+440) = ARR(III+440) + GAUSS3(I) ! Segment 441-660 100 CONTINUE ENDDO C ENDDO C C ------DETERMINE CUTS (I.E. MAX AND MIN VALUES)--------------------- C RMIN = 10000.0 RMAX = 0.0 SCAL = 0.3989422804/SD ! For scaling to given height. DO I = 1, 660 ARR(I) = ARR(I)/SCAL IF (ARR(I).GT.RMAX) RMAX = ARR(I) IF (ARR(I).LT.RMIN) RMIN = ARR(I) ENDDO C C RETURN C END C----------------------------------------------------------------------- SUBROUTINE AG(X1,X2,M,SD,XINTOT) ! Driver for the integr. routine C IMPLICIT NONE C REAL X1,X2,M,SD,INTG1,INTG2 C INTEGER II REAL XINTOT,X,XL,AGAUSS,XU,XINT C XINTOT = 0.0 X = (X2-X1)/100.0 ! Intervals in which integr. will take place. DO II = 1, 100 XL = X1 + FLOAT(II-1)*X XU = X1 + FLOAT(II)*X INTG1 = 100.0*AGAUSS(XL,M,SD) INTG2 = 100.0*AGAUSS(XU,M,SD) IF (ABS(ABS(XL-M)-ABS(XU-M)).LE.1.0E-5) THEN C Here, XL and XU are equally distanced from the mean. XINT = INTG1 ELSE IF (XL.LE.M.AND.XU.GT.M) THEN C Here, XL is on one side, and XU on the other side of the mean. XINT = 0.5*ABS(INTG1-INTG2) + MIN(INTG2,INTG1) ELSE XINT = 0.5*ABS(INTG1-INTG2) ENDIF XINTOT = XINTOT + XINT ENDDO RETURN END C------------------------------------------------------------------------- FUNCTION AGAUSS(X,AVERAG,SIGMA) ! (FROM BEVINGTON, P. 48) C IMPLICIT NONE C DOUBLE PRECISION Z,Y2,TERM,SUM,DENOM,SQR2 C REAL X,AVERAG,SIGMA,AGAUSS C SQR2 = DSQRT(2.D0)/2.0 Z = ABS(X-AVERAG)/SIGMA AGAUSS = 0.0 IF (Z.GT.5.0) THEN ! Precise enough ? AGAUSS = 1.0 RETURN ENDIF IF (Z) 42,42,21 21 TERM = SQR2 * Z SUM = TERM Y2 = (Z**2)/2.0 DENOM = 1.0 31 DENOM = DENOM + 2.0 TERM = TERM*(Y2*2.0/DENOM) SUM = SUM + TERM IF (TERM/SUM -1.0E-10) 41,41,31 41 AGAUSS = 1.128379167*SUM*DEXP(-Y2) 42 RETURN END