C @(#)sgepi.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:21 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 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence 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 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C PROGRAM SGEPI VERSION 1.0 SEP 8, 1988 C P.GROSBOL ESO - GARCHING C C.PURPOSE: C PROGRAM ESTIMATES THE POSITION AND INCLINATION ANGLE OF A GALAXY. C THIS IS DONE BY USE THE VARIATION OF THE 2ND FOURIER HARMONIC. C C.INPUT/OUTPUT: C THE PROGRAM RUNS UNDER MIDAS AND USES THE FOLLOWING KEYWORDS: C IN_A : NAME OF IMAGE FILE (I) C INPUTR(1-5): REAL INPUT PARAMETERS DEFINING RADIAL RANGE (I) C 1:INNER RADIUS, 2:OUTER RADIUS, 3:RADIAL STEP C 4:X CENTER POSITION, 5:Y CENTER POSITION, C OUTPUTR(1-4): REAL OUTPUT PARAMETERS GIVING (O) C 1:X CENTER POSITION, 2:Y CENTER POSITION, C 3:POSITON ANGLE (DEG), 4:INCLINATION ANGLE (DEG) C C.VERSION: C Included in the ESO-MIDAS version 910208 RHW C C 010201 last modif C C--------------------------------------------------------------------- PROGRAM SGEPI C IMPLICIT NONE INTEGER MR INTEGER MF INTEGER MFD INTEGER MIS PARAMETER (MR=100) PARAMETER (MF=6) PARAMETER (MFD=MF+1) PARAMETER (MIS=15) C INTEGER MADRID(1) INTEGER IERR, IAV INTEGER NA INTEGER MC INTEGER I INTEGER NR INTEGER NPIX(3),IMNO,KUNIT(4),INULL INTEGER*8 PNTR REAL PI, PA, PM REAL R, RI, RO, RM, RP REAL DR, DTR, D24, DI REAL XC, YC, XM REAL AN REAL SM, SA, SB REAL WM REAL AIR REAL A, AI REAL B REAL XS, YS, SX, SY C REAL AMP(MFD),PHA(MFD),RANGE(3),POSI(4), 2 AMP2(MR),PHA2(MR),AMP4(MR),PHA4(MR), 3 W(MR),XV(MR),YV2(MR),YV4(MR) DOUBLE PRECISION START(3), STEP(3) C CHARACTER*60 INAME CHARACTER*64 UNIT CHARACTER*72 IDEN CHARACTER*80 MSG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C INITIATE MIDAS C CALL STSPRO('SGEPI') C C READ PARAMETERS C CALL STKRDC('IN_A',1,1,60,IAV,INAME,KUNIT,INULL,IERR) CALL STKRDR('INPUTR',1,3,IAV,RANGE,KUNIT,INULL,IERR) CALL STKRDR('INPUTR',4,2,IAV,POSI,KUNIT,INULL,IERR) C C OPEN IMAGE FRAME C CALL STIGET(INAME,D_R4_FORMAT,F_I_MODE,1,3,NA,NPIX, , START,STEP,IDEN,UNIT,PNTR,IMNO,IERR) IF (NA.NE.2 .OR. IERR.NE.0 .OR. . ABS(STEP(1)/STEP(2)-1.0).GT.1E-5) THEN CALL STTPUT('*** FATAL: WRONG IMAGE FORMAT',IERR) CALL STSEPI ENDIF C C CHECK PARAMETERS AND CONVERT TO PIXEL VALUES C PI = 4.0 * ATAN(1.0) DTR = PI / 180.0 XS = START(1) YS = START(2) SX = STEP(1) SY = STEP(2) RI = MAX(0.0, ABS(RANGE(1))) RO = MAX(RI, ABS(RANGE(2))) RM = MIN(ABS(POSI(1)-XS), ABS(POSI(2)-YS)) RM = MIN(RM, ABS(POSI(1)-XS-NPIX(1)*SX)) RM = MIN(RM, ABS(POSI(2)-YS-NPIX(2)*SY)) RO = MIN(RM-2.0*ABS(SX), RO) DR = ABS(RANGE(3)) XC = (POSI(1)-XS)/SX + 1.0 YC = (POSI(2)-YS)/SY + 1.0 NR = MIN(MR, INT((RO-RI)/DR)) C C GET FOURIER COEFFICIENTS FROM EACH RADIUS C R = RI DO 100, I = 1,NR RP = R / STEP(1) CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC,RP,0.0,0.0, , MFD,MC,AMP,PHA,AN) AMP2(I) = 0.0 PHA2(I) = 0.0 AMP4(I) = 0.0 PHA4(I) = 0.0 W(I) = 0.0 IF (MC.GT.0) THEN AMP2(I) = AMP(3) / AMP(1) PHA2(I) = PHA(3) / (2.0*DTR) AMP4(I) = AMP(5) / AMP(1) PHA4(I) = PHA(5) / (4.0*DTR) D24 = COS( 4.0*DTR*(PHA2(I)-PHA4(I)) ) W(I) = MAX( 0.0 , D24*MIN( 10.0,AMP2(3)/AN) ) ENDIF WRITE(MSG,605) R,AMP2(I),PHA2(I),PHA4(I),D24,AN,W(I) 605 FORMAT(' R,A,P(2,4,6),W :',F6.1,F7.3,2F6.1,2F7.3,F7.3) CALL STTPUT(MSG,IERR) R = R + DR 100 CONTINUE C C COMPUTE ESTIMATE OF POSITION ANGLE C PM = 0.0 RM = 0.0 SM = 0.0 WM = 0.0 DO 200, I = 1,NR RM = RM + W(I)*I XM = XM + W(I)*PHA2(I) WM = WM + W(I) 200 CONTINUE PA = XM / WM RM = (RM/WM - 1.0)*DR + RI WRITE(MSG,620) RM,PA 620 FORMAT('RADIUS, P.A. : ',2F8.1) CALL STTPUT(MSG,IERR) C C ESTIMATE INCLINATION ANGLE WITH THIS POSITION ANGLE C RP = RM / STEP(1) AI = 0.0 DI = 5.0 PA = DTR * PA DO 300, I = 1,MIS AIR = DTR * AI CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC,RP,PA,AIR, , MFD,MC,AMP,PHA,AN) W(I) = 0.0 IF (MC.GT.0) THEN XV(I) = COS( AIR ) YV2(I) = AMP(3) * COS( PHA(3) ) / AMP(1) YV4(I) = AMP(5) * COS( PHA(5) ) / AMP(1) W(I) = 2.0*COS(ATAN(4.0*YV2(I))) ENDIF WRITE(MSG,630) AI,XV(I),YV2(I),YV4(I),W(I) 630 FORMAT(' R,A2,P2,W :',F6.1,F8.3,2F8.3,F8.3) CALL STTPUT(MSG,IERR) AI = AI + DI 300 CONTINUE C C FIT LINE TO GET ESTIMATE FOR INCLINATION ANGLE C CALL LFIT(XV,YV2,W,MIS,1,A,SA,B,SB,R) AI = ACOS( -A/B ) / DTR PA = PA / DTR WRITE(MSG,640) R,PA,AI 640 FORMAT('RADIUS, PA, I : ',F8.1,2F7.1) CALL STTPUT(MSG,IERR) CALL STKWRR('OUTPUTR',POSI,1,2,KUNIT,IERR) CALL STKWRR('OUTPUTR',PA,3,1,KUNIT,IERR) CALL STKWRR('OUTPUTR',AI,4,1,KUNIT,IERR) C C FINISHED - EXIT C CALL STSEPI END