SUBROUTINE PSF_RAD (R,EN,NEBN,ISTR,FRAD,STATUS)
C
C------------------------------------------------------------------------
C
C.IDENTIFICATION:
C	Subroutine psf_rad
C
C.SOURCE FILE:
C       psf_rad.f	
C
C.LIBRARY:
C       libmecs.a
C
C.AUTHOR:
C       S.Molendi    IFCTR/CNR,  Milano
C
C.VERSIONS:
C	Vers. 1.0     6-7-95		creation date
C       Vers. 1.1    24-7-95            revised by T.Mineo
C       Vers. 2.1    14-5-96            revised by T.Mineo
c       Vers. 2.2    17-5-96            XASized by L.Chiappetti
C
C.PURPOSE:
C       computes the  correction to the total efficiency due to the 
C       selection inside a spot of R dimension
C
C.METHOD
C       This subroutine uses the integral on-axis PSF model 
C       to compute the fraction of counts falling
C       inside a given extraction radius R vs Energy.
C       The shape of the PSF is assumed not changing for off-axis 
C       source.
C       The program normalizes the PSF to the values at 15 mm (ie
C       FRAD(15mm)=1)
C
C.SYNTAX:
C	CALL PSF_RAD (R,EN,NEBN,ISTR,FRAD,STATUS)
C
C.PARAMETERS:
C       inp. R          Real*4    Selection radius (mm)
C       inp. EN(NEBN)   Real*4    energy (keV)
C       inp. NEBN       Integer*4 number of energy bins
C       inp. ISTR       Integer*4 1=ME1, 2=ME2, 3=ME3
C       out. FRAD(NEBN) Real*4    efficiency correcting factor
C       out. STATUS     Integer*4 STATUS=0 ==>no errors 
C
C.REFERENCES:
C       S.Molendi  Sax technical report .......
C
C------------------------------------------------------------------------
C
C...variables definition
C
	INTEGER*4 ISTR,STATUS,NEBN,I,IFLAG_PSF
C
	REAL*4 EN(NEBN),FRAD(NEBN),R
        REAL*4 A_SIGMA(3),B_SIGMA(3),C_SIGMA(3),D_SIGMA(3),
     .         E_SIGMA(3)
        REAL*4 A_RL(3),B_RL(3),C_RL(3),D_RL(3),E_RL(3)
        REAL*4 A_M(3),B_M(3),C_M(3),D_M(3),E_M(3)
        REAL*4 A_RAT(3),B_RAT(3),C_RAT(3),D_RAT(3),E_RAT(3)
	REAL*4 PS(3)
        REAL*4 SIGMA,RL,M,RAT,R_DET,E
        REAL*4 INT_GAUSS,INT_GEN_L,NOR_GAUSS,NOR_GEN_L
C
        CHARACTER BUFFER*80,NAME*80
        INTEGER   ISYS
        LOGICAL   VOSERROR
C
C...R_DET=detector radius in mm
C
	DATA R_DET /15.0/
C
	SAVE IFLAG_PSF
C
C...checks the rsel radius
C
	IF(R.GT.R_DET) THEN
	  WRITE(*,*)'PSF_RAD: spot selection radius > detector radius'
	  STATUS=1
	  RETURN
	END IF
C
C...reads the input parameters only once
C
	IF(IFLAG_PSF.EQ.0) THEN
        WRITE(BUFFER, 120)ISTR
 120    FORMAT('m',I1,'_psf.coeff')
        CALL BUILDPATH(BUFFER,'CALIB',NAME)
        CALL Z_OPEN(12,NAME,'SEQ','OLD',0)
        IF (VOSERROR(STATUS,ISYS))RETURN
C
C...Parametres for the Sigma of the Gaussian Component
C
	  READ(12,*) 
	  READ(12,*) 
          READ(12,*) A_SIGMA(ISTR)
          READ(12,*) B_SIGMA(ISTR)
          READ(12,*) C_SIGMA(ISTR)
          READ(12,*) D_SIGMA(ISTR)
          READ(12,*) E_SIGMA(ISTR)
C
C...Parameters for the radius of the Generalized Lorenzian
C
	  READ(12,*) 
          READ(12,*) A_RL(ISTR)
          READ(12,*) B_RL(ISTR)
          READ(12,*) C_RL(ISTR)
          READ(12,*) D_RL(ISTR)
          READ(12,*) E_RL(ISTR)
C
C...Parameters for the index of the Generalized Lorenzian
C
	  READ(12,*) 
          READ(12,*) A_M(ISTR)
          READ(12,*) B_M(ISTR)
          READ(12,*) C_M(ISTR)
          READ(12,*) D_M(ISTR)
          READ(12,*) E_M(ISTR)
C
C...Ratio of the Gaussian Normalization over the Generalized Lorenzian 
C
	  READ(12,*) 
          READ(12,*) A_RAT(ISTR)
          READ(12,*) B_RAT(ISTR)
          READ(12,*) C_RAT(ISTR)
          READ(12,*) D_RAT(ISTR)
          READ(12,*) E_RAT(ISTR)
C
C...reads the pixel platescale 
C
	  READ(12,*)
          READ(12,*) PS(ISTR)
C
	  CLOSE(12)
	  IFLAG_PSF=1
	END IF
C
	DO 10 I=1,NEBN
	  E=EN(I)
C
C...Computation of the parameters
C
          SIGMA = PS(ISTR)*(C_SIGMA(ISTR)*
     .        EXP(-(E-A_SIGMA(ISTR))/B_SIGMA(ISTR))+
     .        D_SIGMA(ISTR)+E_SIGMA(ISTR)*E)
          RL = PS(ISTR)*(C_RL(ISTR)*EXP(-(E-A_RL(ISTR))/B_RL(ISTR))+
     .         D_RL(ISTR)+E_RL(ISTR)*E)
          M   = C_M(ISTR)*EXP(-(E-A_M(ISTR))/B_M(ISTR))+D_M(ISTR)+
     .        E_M(ISTR)*E
          RAT = C_RAT(ISTR)*EXP(-(E-A_RAT(ISTR))/B_RAT(ISTR))+
     .        D_RAT(ISTR)+E_RAT(ISTR)*E
C
C...Computation of the Integral PSF
C
          INT_GAUSS = RAT*SIGMA*SIGMA*
     .                (1.-EXP(-1.*R*R/2./SIGMA/SIGMA)) 
          INT_GEN_L = RL*RL/2./(M-1.)*
     .                (1.-1./(1.+(R/RL)**2)**(M-1.))
          NOR_GAUSS = RAT*SIGMA*SIGMA*
     .                (1.-EXP(-1.*R_DET*R_DET/2./SIGMA/SIGMA)) 
          NOR_GEN_L = RL*RL/2./(M-1.)*
     .                (1.-1./(1.+(R_DET/RL)**2)**(M-1.))
          FRAD(I)   = (INT_GAUSS+INT_GEN_L)/(NOR_GAUSS+NOR_GEN_L)
C
 10	CONTINUE
C
      RETURN
      END