C @(#)rfotmodel.for	19.1 (ES0-DMD) 02/25/03 13:30:15
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 MODEL
C+++
C.PURPOSE: Interpolation program for undersampling problem in ROMAFOT
C.AUTHOR:  R. Buonnanno
C.REVISED: R.H. Warmels
C.VERSION:  
C---------------
      IMPLICIT         NONE
C
      INTEGER          MTG(55,55)
      DOUBLE PRECISION DX,DY
      DOUBLE PRECISION IN,ER
      DOUBLE PRECISION DX0,DY0
      DOUBLE PRECISION VOLU,VOLG,TNOIS
      INTEGER          IAV, ISTAT
      INTEGER          IDX, IDY
      INTEGER          I, J, II
      INTEGER          INCR
      INTEGER          KUN, KNUL
      INTEGER          KVX, KVY
      INTEGER          KVM1, KVM2, KVM3, KVM
      INTEGER          MC, MR
      INTEGER          NCHAR
      INTEGER          NP3, NPUG, NSA
      INTEGER          NPX, NPY
      INTEGER          NTOT, NPXS, NPYS
C
      REAL             AMD, AML
      REAL             BAF
      REAL             BET, DITE, EFF, FCON
      REAL             COSTMA
      REAL             ERROLD, ERR
      REAL             FX, FY, FLUSSO, FSP
      REAL             H
      REAL             PMA, PPO, PG, PSF
      REAL             QS
      REAL             RINPUT(15)
      REAL             RMA, RON, RNOIS
      REAL             SEEING, SIG, SKYL
      REAL             TESP, TET
      REAL             XB, YB
      REAL             X0, Y0
      REAL             XA, YA
      REAL             VM, VX, VY, VT, VF
      REAL             VFCG, VFCGS
      REAL             GALE2
C
      CHARACTER        FRE*1,STRING*80,FILE*60
C
      DATA             PG/3.14159/
      DATA             SKYL/3.18E-06/
C
  999 FORMAT('*** INFO: Ratio FWHM over Pixel size =', E11.6)
  998 FORMAT(2X,I3,2X,I3,2X,A1)
  997 FORMAT(2X,F5.3,3(2X,F6.2),2X,F9.4)
  996 FORMAT(2(2X,F3.1),2X,F6.2)
  995 FORMAT(25I3)
C
C *** See in MIDAS a alive
      CALL STSPRO('MODEL')
C
C *** the default setting:
C                           DEFAULT:
C  SKYL=3.18E-06            photon surface brightness of the sky in 
C                           photons cm**-2 s**-1 A**-1 arcsec**-2
C                           corresponding to V=21.3 at 5500 A
C  SEEING1 = 1              seeing (in arcsec)
C  BET1    = 3              beta
C  DITE1   = 2.2            diameter of the telescope (metres)
C  BAF1    = 800.           bandwidth
C  EFF1    = .3             total efficiency 
C  FCON1   = 10             photons-to-ADU conversion factor
C  TESP1   = 1000.          exposure time
C  TET1    = .35            pixel size (in arcsec) 
C  RON     = 35             read out noise of the detector  
C  AMD1    = 23.            limiting magnitude (faint end)
C  AML1    = 18.            limiting magnitude (bright end)
C  PMA1    = .5             sampling in magnitude
C  PPO     = .1             positional accuracy??
C
C *** get the input parameters
      CALL STKRDR('INPUTR',1,13,IAV,RINPUT,KUN,KNUL,ISTAT)
      SEEING = RINPUT(1)                                            ! seeing
      BET    = RINPUT(2)
      DITE   = RINPUT(3)*100.
      BAF    = RINPUT(4)
      EFF    = RINPUT(5)
      FCON   = RINPUT(6)
      TESP   = RINPUT(7)
      TET    = RINPUT(8)
      RON    = RINPUT(9)
      AMD    = RINPUT(10)
      AML    = RINPUT(11)
      PMA    = RINPUT(12)
      QS     = RINPUT(13)
C
      FSP    = SEEING/TET
      SIG    = (SEEING/TET)/(2*SQRT(2**(1./BET)-1))
      PPO    = .1
      COSTMA = -2.5*ALOG10(FCON/(DITE**2*BAF*TESP*EFF*PG/4.))+7.5
      MR     = 7+INT(FSP/1.4+.7)*2
      IF (MR.GT.15) THEN
         MR = 15
      ENDIF
      MC     = MR
      IDX    = 0
      IDY    = 0
      X0     = MC/2
      Y0     = MR/2
      FRE    = 'N'
      XA     = INT(X0)-MC/2+IDX
      YA     = INT(Y0)-MR/2+IDY
      XB     = XA+MC
      YB     = YA+MR
C
C ***
      WRITE(STRING,999) FSP
      CALL STTPUT(STRING,ISTAT)
      CALL STKRDC('INPUTC',1,1,60,IAV,FILE,KUN,KNUL,ISTAT)
      NCHAR = INDEX(FILE,' ')-1
      FILE = FILE(1:NCHAR)//'.dat'
      OPEN(UNIT=1,STATUS='NEW',ACCESS='SEQUENTIAL',FORM='FORMATTED',
     2     FILE=FILE)
C
C *** write the output table
      WRITE(1,998) MR,MC,FRE
      WRITE(1,997) PPO,AMD,AML,PMA,COSTMA
C
      INCR = -INT(PPO*100+0.5)
      DO 1000 KVX = 50,0,INCR
         VX   = FLOAT(KVX)/100.
         X0   = MC/2+VX
         INCR = -INT(PPO*100+0.5)
         DO 1100 KVY = 50,KVX,INCR
            VY = FLOAT(KVY)/100.
            Y0 = MR/2+VY
            KVM1 = INT(100*AML+0.5)
            KVM2 = INT(100*AMD+0.5)
            KVM3 = INT(100*PMA+0.5)
            DO 1200 KVM = KVM1,KVM2,KVM3
               VM     =  FLOAT(KVM)/100.
               RMA    = VM
               FLUSSO = 10.**(-(RMA-7.5)/2.5)
               H      = (BET-1)*FLUSSO*DITE**2*BAF*TESP*EFF/
     2                  (4*SIG**2*FCON)
               FX     = X0-INT(X0)
               FY     = Y0-INT(Y0)
               WRITE(1,996) FX,FY,RMA
C
C *** calcuate matrix
               NP3   = 0
C              NPUG  = NPU
               VOLU  = 0.
               VOLG  = 0.
               TNOIS = 0.
               DO 1300 I=1,MR
                  YA=INT(Y0)-MR/2+I-1+IDY
                  DO 1400 J=1,MC
                     XA=INT(X0)-MC/2+J-1+IDX
C
C *** calcualte int vero
                     IF (BET.EQ.1.5) THEN
                        DX  = XA+1-X0
                        DY  = YA+1-Y0
                        DX0 = DX-1.
                        DY0 = DY-1.
                        ER  = SIG
                        VT  = H*(IN(DX,DY,ER)-IN(DX,DY0,ER)-
     2                          IN(DX0,DY,ER)+IN(DX0,DY0,ER))*ER**2
                     ELSE
                        VT  = GALE2(4,4,H,SIG,BET,XA,YA,X0,Y0)
                     END IF
                     FLUSSO = VT*FCON/(EFF*DITE**2*BAF*TESP*PG/4.)
                     RNOIS  = SQRT((FLUSSO+SKYL*TET**2.)*EFF*TESP*
     2                        BAF*DITE**2.*PG/4+RON**2)
C
                     TNOIS  = TNOIS+RNOIS**2
                     VOLU   = VOLU+VT
                     VF     = PSF(H,XA+.5,YA+.5,X0,Y0,SIG,BET)
C
                     IF (FRE.EQ.'N') THEN
                        VFCG = VF
                        ERR  = RNOIS*QS/FCON
                        NPUG = 1
 1401                   IF (ABS(VT-VFCG).LE.ERR.OR.NPUG.GT.5) GO TO 1402
                           NPUG = NPUG+1
                           VFCG = GALE2(NPUG,NPUG,H,SIG,BET,XA,YA,X0,Y0)
                           GO TO 1401
 1402                   CONTINUE
                        VOLG     = VOLG+VFCG 
                        MTG(I,J) = NPUG**2
                        NP3      = NP3+NPUG**2
                     ELSE
                        ERR      = RNOIS*QS/FCON
                        NSA      = 1000
                        ERROLD   = 10.E15
                        DO 1500 NPX=1,4
                           DO 1600 NPY=1,4
                              VFCG=GALE2(NPX,NPY,H,SIG,BET,XA,YA,X0,Y0)
                              IF (ABS(VT-VFCG).LT.ERR) THEN
                                 NTOT=NPX*NPY
                                 IF (NTOT.LE.NSA) THEN
                                    IF (ERR.LT.ERROLD. OR. 
     2                                  NTOT.NE.NSA) THEN
                                        NSA    = NTOT
                                        NPXS   = NPX
                                        NPYS   = NPY
                                        VFCGS  = VFCG
                                        ERROLD = ERR
                                     END IF
                                  END IF
                              END IF
 1600                      CONTINUE
 1500                   CONTINUE
C
                        VFCG     = VFCGS
                        VOLG     = VOLG+VFCG
                        MTG(I,J) = 10*NPXS+NPYS
                        NP3      = NP3+NSA
                     END IF 
 1400             CONTINUE
 1300          CONTINUE
C 
               DO 1201 II=1,MR
                  WRITE(1,995) (MTG(II,J),J=1,MC)
 1201          CONTINUE
 1200       CONTINUE
 1100    CONTINUE
 1000 CONTINUE
C
      CLOSE(UNIT=1)
      CALL STSEPI
      END


      REAL FUNCTION PSF(H,X,Y,X0,Y0,SIG,BET)
C
      REAL  H
      REAL  X
      REAL  Y
      REAL  X0
      REAL  Y0
      REAL  SIG
      REAL  BET
C
      PSF = H * (1+((X-X0)**2 + (Y-Y0)**2)/SIG**2)**(-BET)
C
      RETURN
      END


      DOUBLE PRECISION FUNCTION IN(DX,DY,ER)
C
      DOUBLE PRECISION DX
      DOUBLE PRECISION DY
      DOUBLE PRECISION ER
C
      IN = DATAN(DX*DY/(ER*DSQRT(DX**2+DY**2+ER**2)))
C
      RETURN
      END


