C @(#)surface.for	19.1 (ES0-DMD) 02/25/03 14:25:41
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 SURFACE
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT        (c) 1988 European Southern Observatory
C.IDENTIFICATION   SURFACE
C.LANGUAGE         ESO-FOR
C.AUTHOR           Preben Grosbol (ESO-IPG)
C.KEYWORDS         Artificial objects, simulation
C.ENVIRONMENT      MIDAS
C.COMMENT          MIDAS command similar to the IHAP SURFACE command
C.PURPOSE          The command add an artificial object or mathematical
C                  surface to an existing frame
C.VERSION          1.0    1989-Feb-13 : Creation      (PJG)
C---------------------------------------------------------------------
C
      IMPLICIT NONE
      INTEGER MTROW
      PARAMETER  (MTROW=512)
C
      CHARACTER*2       METH
      CHARACTER*40      INAME,TNAME
      INTEGER           IAV,IERR,NA,ILAB
      INTEGER           NPIX(2),KUNIT(4),INULL
      INTEGER           IXF,IXL,IYF,IYL,I,I1,I2,NX,NY
      INTEGER           IDF,IMEM,IDX,IOFF,IA,IB,IX,IY,IVA
      REAL              CI,CP,SP,X,Y,XS,YS,DX,DY,XC,YC,XT,YT
      REAL              RL,R,SC,AMP,EL,ELIM,RQCNST,DTR
      REAL              FX,FY,FXX,FYY,FXY,FXXX,FYYY,FXXY,FYYX
      REAL              GA2RAN
      REAL              POSI(2),PROJ(2),PARM(10)
      REAL              VMADD(1)
C     REAL              SAMP(MTROW),SPHA(MTROW)
      DOUBLE PRECISION  START(2),STEP(2)
C
      COMMON /VMR/VMADD
C
      EQUIVALENCE (NX,NPIX(1)),(NY,NPIX(2))
      EQUIVALENCE (XS,START(1)),(YS,START(2))
      EQUIVALENCE (DX,STEP(1)),(DY,STEP(2))
      EQUIVALENCE (XC,POSI(1)),(YC,POSI(2)),(AMP,PARM(1))
C
      INCLUDE 'MID_REL_INCL:MIDAS_DEF.INC'
C
C     INITIATE MIDAS
C
      CALL STSPRO('SURFACE ',IERR)
C
C     INITIATE VARIABLES
C
      ELIM = -20.0
      RQCNST = -3.33
C
C     READ FRAME NAME AND OPEN IT
C
      CALL STKRDC('IN_A ',1,1,40,IAV,INAME,KUNIT,INULL,IERR)
      CALL STFOPN(INAME,D_R4_FORMAT,0,1,IDF,IERR)
C
C     READ FRAME DESCRIPTORS
C
      CALL STDRDI(IDF,'NAXIS ',1,1,IAV,NA,KUNIT,INULL,IERR)
      IF (NA.GT.2 .OR. IERR.NE.0) THEN
         CALL STTPUT('WRONG IMAGE DIMENSION >2 ',IERR)
         GOTO 90000
      ENDIF
      CALL STDRDI(IDF,'NPIX ',1,2,IAV,NPIX,KUNIT,INULL,IERR)
      CALL STDRDD(IDF,'START ',1,2,IAV,START,KUNIT,INULL,IERR)
      CALL STDRDD(IDF,'STEP ',1,2,IAV,STEP,KUNIT,INULL,IERR)
C
C     READ PARAMETERS 
C
      CALL STKRDR('POSITION ',1,2,IAV,POSI,KUNIT,INULL,IERR)
      CALL STKRDR('PARM ',1,10,IAV,PARM,KUNIT,INULL,IERR)
C
C     INITIATE VARIABLES
C
      IF (NA.LE.1) THEN
         NY = 1
         YS = 0.0
         DY = 1.0
         YC = 0.0
      ENDIF
      IYF = 1
      IYL = NY
      IXF = 1
      IXL = NX
      EL = ALOG(10.0) * ELIM
C
C     READ AND CHECK METHOD
C
      CALL STKRDC('METH ',1,1,2,IAV,METH,KUNIT,INULL,IERR)
      IF (METH.EQ.'PO') THEN
         ASSIGN 1000 TO ILAB
         FX = PARM(2)
         FY = PARM(3)
         FXX = PARM(4)
         FYY = PARM(5)
         FXY = PARM(6)
         FXXX = PARM(7)
         FYYY = PARM(8)
         FXXY = PARM(9)
         FYYX = PARM(10)
      ELSEIF (METH.EQ.'GA') THEN
         ASSIGN 2000 TO ILAB
         IF (PARM(2).LE.0.0) PARM(2) = 1.0E-10
         SC = -2.772589/(PARM(2)*PARM(2))
         RL = 0.7 * SQRT(EL/SC)
         IXF = MAX(1,MIN(NX,INT((XC-RL-XS)/DX)+1))
         IXL = MAX(1,MIN(NX,INT((XC+RL-XS)/DX)+1))
         IYF = MAX(1,MIN(NY,INT((YC-RL-YS)/DY)+1))
         IYL = MAX(1,MIN(NY,INT((YC+RL-YS)/DY)+1))
      ELSEIF (METH.EQ.'DI') THEN
         ASSIGN 3000 TO ILAB
         IF (PARM(2).LE.0.0) PARM(2) = 1.0E-10
         SC = -1.0/PARM(2)
         RL = 0.7 * EL/SC
         IXF = MAX(1,MIN(NX,INT((XC-RL-XS)/DX)+1))
         IXL = MAX(1,MIN(NX,INT((XC+RL-XS)/DX)+1))
         IYF = MAX(1,MIN(NY,INT((YC-RL-YS)/DY)+1))
         IYL = MAX(1,MIN(NY,INT((YC+RL-YS)/DY)+1))
      ELSEIF (METH.EQ.'EL') THEN
         ASSIGN 4000 TO ILAB
         IF (PARM(2).LE.0.0) PARM(2) = 1.0E-10
         SC = 1.0/PARM(2)
         RL = 0.7 * ((ELIM/RQCNST)**4)/SC
         IXF = MAX(1,MIN(NX,INT((XC-RL-XS)/DX)+1))
         IXL = MAX(1,MIN(NX,INT((XC+RL-XS)/DX)+1))
         IYF = MAX(1,MIN(NY,INT((YC-RL-YS)/DY)+1))
         IYL = MAX(1,MIN(NY,INT((YC+RL-YS)/DY)+1))
      ELSEIF (METH.EQ.'SP') THEN
         ASSIGN 5000 TO ILAB
         CALL STKRDC('IN_B ',1,1,40,IAV,TNAME,KUNIT,INULL,IERR)
      ELSEIF (METH.EQ.'NO') THEN
         ASSIGN 6000 TO ILAB
         IF (0.0.GT.PARM(2)) THEN
            IA = PARM(2)
            IB = PARM(3)
            IX = PARM(4)
         ELSE
            IA = 173
            IB = 3297
            IX = 11321
         ENDIF
         CALL USEED(IA,IB,IX)
      ELSE
         CALL STTPUT('METHOD NOT SUPPORTED ',IERR)
         GOTO 90000
      ENDIF
C
C     READ AND INITIATE PROJECTION 
C
      CALL STKRDR('PROJECT ',1,2,IAV,PROJ,KUNIT,INULL,IERR)
      DTR = ATAN(1.0D0)/45.0D0
      CI = COS(DTR*PROJ(1))
      CP = COS(DTR*PROJ(2))
      SP = SIN(DTR*PROJ(2))
C
C     GET MEMORY AND INITIATE VARIABLES
C
      CALL STFCRE('DUMMY ',D_R4_FORMAT,F_X_MODE,1,NX,IMEM,IERR)
      CALL STFMAP(IMEM,F_X_MODE,1,NX,IVA,IDX,IERR)
      IOFF = NX*(IYF-1) + 1
C
C     GO THROUGH FRAME LINE BY LINE
C
      Y = YS - YC + (IYF-1)*DY
      DO 100, IY = IYF,IYL
         CALL STFGET(IDF,IOFF,NX,IVA,VMADD(IDX),IERR)
         I1 = IDX + IXF - 1
         I2 = IDX + IXL - 1
         X = XS - XC + (IXF-1)*DX
         DO 300, I = I1,I2
            GOTO ILAB
C
C     POLYNOMIUM
C
 1000       CONTINUE
            VMADD(I) = VMADD(I) + AMP + X*(FX+X*(FXX+X*FXXX+Y*FXXY)) +
     +                 Y*(FY+X*FXY+Y*(FYY+Y*FYYY+X*FYYX))
            GOTO 200
C
C     GAUSSIAN
C
 2000       CONTINUE
            R = X*X + Y*Y
            VMADD(I) = VMADD(I) + AMP*EXP(SC*R)
            GOTO 200
C
C     EXPONENIAL DISK
C
 3000       CONTINUE
            XT = SP*X - CP*Y
            YT = CI * (CP*X + SP*Y)
            R = SQRT(XT*XT + YT*YT)
            VMADD(I) = VMADD(I) + AMP*EXP(R*SC)
            GOTO 200
C
C     ELLIPTICAL GALAXY : R**1/4 LAW
C
 4000       CONTINUE
            XT = SP*X - CP*Y
            YT = CI * (CP*X + SP*Y)
            R = SQRT(XT*XT + YT*YT)
            VMADD(I) = VMADD(I) + AMP*10.0**(RQCNST*(R*SC)**(0.25))
            GOTO 200
C
C     SPIRAL PATTERN
C
 5000       CONTINUE
            XT = SP*X - CP*Y
            YT = CI * (CP*X + SP*Y)
            R = SQRT(XT*XT + YT*YT)
            VMADD(I) = VMADD(I) + AMP
            GOTO 200
C
C     GAUSSIAN NOISE
C
 6000       CONTINUE
            VMADD(I) = VMADD(I) + AMP*GA2RAN()
C
  200       CONTINUE
            X = X + DX
  300    CONTINUE
         Y = Y + DY
         CALL STFPUT(IDF,IOFF,NX,VMADD(IDX),IERR)
         IOFF = IOFF + NX
  100 CONTINUE
C
C     FINISHED - EXIT
C
90000 CALL STFCLO(IDF,IERR)
      CALL STFCLO(IMEM,IERR)
      CALL STSEPI(IERR)
      END
