C @(#)surface.for 17.1.1.1 (ES0-DMD) 01/25/02 17:54:07 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