C @(#)spec3g.for	19.1 (ES0-DMD) 02/25/03 13:24:48
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       
