c*********************************************************************
c     calculates Compton reflected spectrum 
c     from the angular dependent incoming radiation 
c     using Green's matrix from Poutanen, Nagendra & Svensson 1996 
c can be used only for NEUTRAL reflection 
c*********************************************************************
      subroutine PNSrefl(compsp,xr,frefl,angout,
     + par_refl,nrefl,jref,kko,ncout,ioniz,irefl)
c     XX(II) - energies where comptonized spectrum is given 
c     compsp(kk) - angle dependent comptonized spectrum (intensities)
c     UANG(NC)  -angles where comptonized spectrum is given 
c     xr(jref) - energies where frefl should be computed
c     frefl(kko)  - reflected spectrum (intensities) at different angles
c     angout(ncout) - angles where  frefl should be computed
      PARAMETER(nem=5000) 
      PARAMETER(MAXFRE=71,MAXANG=3)
      PARAMETER(II=MAXFRE,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK) 
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NRG=7,NFEAB=3,NCA=3)
      PARAMETER(NCSM=50,KKC=NCA*NCSM,LCSM=ND*KKC)
      PARAMETER(NCRG=NRG*NCSM,KCRG=NCA*NCRG)
      PARAMETER(IPE=NCSM+NCSM*(NCSM-1)/2,KPE=NCA*IPE)

      implicit real*8 (a-h,o-z)
      integer jref,kko,ncout,ioniz,irefl 
      integer i,k,nrefl
      real*8 compsp(kk),specold(kk) 
      real*8 xr(jref),frefl(kko),angout(ncout)
      real*4 cmin,par_refl(nrefl),ab_,ab0,ab0old  
      real*8 dlamb(nem)
      COMMON/QQCM2/A(II),UANG(NC),AANG(NC),AC(LL),AINT(LL),CINT(LL)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/CIN/UANIN(NCA),AANIN(NCA)
      COMMON/RCSMA/RGCSM(KPE,NCA,NRG,NFEAB),XCSM(NCSM,NRG)
      COMMON/REFEAB/FEABREF(NFEAB)
      DIMENSION FCIN(KKC),FCOUT(KKC)

      DIMENSION FXREF(KCRG),DLXCSM(NCSM),DCSXX(NCSM)
      DIMENSION DLXCRG(NCRG),XCRG(NCRG),FCREF(KCRG)
      DIMENSION FINT(NC),FOUT(NC) 

      logical flagpns,feflag,flagsp 
      save flagpns,feflag,flagsp 
      save ab0old,FCREF,DLXCRG,XCRG,ICO,specold
      data cmin/1d-30/,ab0old/999./
      data IRGMIN/2/ 
      data flagpns/.true./,feflag/.true./,flagsp/.true./  
c Fe abundance 
       ab0=par_refl(1) 
       FEABLG=dble(ab0) 
c metal abundance 
       ab_=par_refl(2) 

       if(ab0.ne.ab0old) then 
       ab0old=ab0 
       feflag=.true. 
       endif 

       do k=1,kk 
       if(compsp(k).ne.specold(k)) then 
       flagsp=.true.
       specold(k)=compsp(k)
       endif 
       enddo 

      if(flagpns) then 
c matrix for Compton reflection by Poutanen, Nagendra & Svensson
      CALL QDRGSDO(UANIN,AANIN,NCA)
      DO  L=1,NCA 
         UANIN(L)=0.5D0*(1D0+UANIN(L))
         AANIN(L)=AANIN(L)*5D-1
      ENDDO 
      CALL REFLMATR
      flagpns=.false.
      endif
      CALL DLOCATE(FEABREF,NFEAB,FEABLG,IFE)

c computation of lambda at xr 
      call DLAMBDA(ab_,ab0,xr,DLAMB,jref) 

c compute some of that only  if FEAB is changed or compsp changed 
      if(feflag.or.flagsp) then 
C calculation of reflected spectrum at XCSM
      DO 100 IRG=IRGMIN,NRG

      DO 110 IANG=1,NC
      DO 110 I=1,NCSM
      KI=IANG+(I-1)*NC
      CALL QHUNT(XX,II,XCSM(I,IRG),ICS)
      IF(ICS.lt.1.or.ICS.ge.II) THEN
      FCIN(KI)=CMIN
      ELSE
c interpolation of Comptonized radiation (down) to XCSM
      CALL DLINEAR(XLOG(ICS),DLOG(COMPSP(IANG+(ICS-1)*NC)),
     *  XLOG(ICS+1),DLOG(COMPSP(IANG+ICS*NC)),
     *  DLOG10(XCSM(I,IRG)),YYY) 
      FCIN(KI)=DEXP(YYY)
      ENDIF
  110 CONTINUE


      CALL REFLCOMP(FCIN,FCOUT,IFE,IRG)

c interpolation of FCOUT to XCSM(1,IRG)
c and renormalization
      IF(IRG.gt.IRGMIN) THEN
      DO 114 IANG=1,NC
      DO 116 I=1,NCSM
      K=IANG+(I-1)*NC
      KRG=K+(IRG-2)*KKC
      DLXCSM(I)=DLOG(XCSM(I,IRG-1))
      DCSXX(I)=DLOG(FXREF(KRG)+CMIN)
  116 CONTINUE 
      XCFLG=DLOG(XCSM(1,IRG)) 
      CALL DLOCATE(DLXCSM,NCSM,XCFLG,ICS)
      CALL DLINEAR(DLXCSM(ICS),DCSXX(ICS), 
     *  DLXCSM(ICS+1),DCSXX(ICS+1),
     *  XCFLG,YYY) 
      FOUT(IANG)=DEXP(YYY)/FCOUT(IANG)
  114 CONTINUE
      ELSE
      DO 115 IANG=1,NC
      FOUT(IANG)=1D0
  115 CONTINUE
      ENDIF

      DO 118 IANG=1,NC
      DO 118 I=1,NCSM
      K=IANG+(I-1)*NC
      KRG=K+(IRG-1)*KKC
      FXREF(KRG)=FCOUT(K)*FOUT(IANG)
  118 CONTINUE

  100 CONTINUE

c calculation of the reflection on a long grid
      ICO=1
      DO IRG=IRGMIN,NRG-1
       I=1
       DO WHILE(XCSM(I,IRG).lt.XCSM(1,IRG+1).and.I.le.NCSM)
       XCRG(ICO)=XCSM(I,IRG)
       DLXCRG(ICO)=DLOG10(XCSM(I,IRG))
        DO J=1,NC
        K=J+(I-1)*NC
        KCO=J+(ICO-1)*NC
        KRG=K+(IRG-1)*KKC
        FCREF(KCO)=FXREF(KRG)
        ENDDO
       I=I+1
       ICO=ICO+1
       ENDDO
      ENDDO
      DO I=1,NCSM
       XCRG(ICO)=XCSM(I,NRG)
       DLXCRG(ICO)=DLOG10(XCSM(I,NRG))
        DO J=1,NC
        K=J+(I-1)*NC
        KCO=J+(ICO-1)*NC
        KRG=K+(NRG-1)*KKC
        FCREF(KCO)=FXREF(KRG)
        ENDDO
       ICO=ICO+1
      ENDDO

      ICO=ICO-1

      feflag=.false. 
      flagsp=.false. 
      endif 
c endif for feflag & flagsp 

c calculation of reflection at out grid
      XCSMIN=XCSM(1,IRGMIN)
      CALL DLOCATE(XR,JREF,XCSMIN,ICS)

c Compton reflection 
      DO 128 I=ICS+1,JREF 
      CALL QHUNT(XCRG,ICO,XR(I),IRE)
      DO 127 IANG=1,NC 
      IF(IRE.lt.1.or.IRE.ge.ICO) THEN
      FINT(IANG)=CMIN
      ELSE
      CALL DLINEAR(DLXCRG(IRE),DLOG(FCREF(IANG+(IRE-1)*NC)+CMIN),
     *  DLXCRG(IRE+1),DLOG(FCREF(IANG+IRE*NC)+CMIN),
     *  DLOG10(XR(I)),YYY)
      FINT(IANG)=DEXP(YYY)
      ENDIF
  127 CONTINUE
c interpolation to the needed angles 
      DO 128 J=1,NCOUT 
      K=J+(I-1)*NCOUT 
      CALL POLINTQ(UANG,FINT,NC,ANGOUT(J),FLY,DELY)
      FREFL(K)=FLY 
  128 CONTINUE

c  Rayleigh Reflection
c here it is dangerous if xr(ics+1) is much bigger than 5 keV 
      DO 136 I=1,ICS+1
      CALL QHUNT(XX,II,XR(I),IRE)
      DO 138 IANG=1,NC
      IF(IRE.lt.1.or.IRE.ge.II) THEN
      FINT(IANG)=CMIN
      ELSE
      CALL DLINEAR(XLOG(IRE),DLOG(COMPSP(IANG+(IRE-1)*NC)),
     *  XLOG(IRE+1),DLOG(COMPSP(IANG+IRE*NC)),
     *  DLOG10(XR(I)),YYY)
      FINT(IANG)=DEXP(YYY)
      ENDIF 
  138 CONTINUE
      CALL RAYREFL(FINT,DLAMB(I),FOUT)


c interpolation to the needed angles
      DO 136 J=1,NCOUT
      K=J+(I-1)*NCOUT
      CALL POLINTQ(UANG,FOUT,NC,ANGOUT(J),FLY,DELY)
      if(I.le.ICS) then
      FREFL(K)=FLY
      else
c use value from Compton reflection FREFL(J+ICS*NC) 
      FINT(J)=FLY/FREFL(K)
      endif
  136 CONTINUE

c renormalization of the Compton part
      do I=ICS+1,JREF 
       do J=1,NCOUT 
       K=J+(I-1)*NCOUT 
       FREFL(K)=FREFL(K)*FINT(J)
       enddo
      enddo

      return
      end


c**************************************************
      SUBROUTINE DLAMBDA(ab_,ab0,XX,DLAMB,II) 
      integer ix,ii,nem 
      parameter(nem=5000) 
      real*8 XX(II),DLAMB(II),eps 
      real*8 X,CS,abscon,dla    
      real*4 xr(nem),xnor,sigmaneu,ab_,ab0
      data eps/1d-10/

      DO IX=1,II
      xr(ix)=sngl(XX(IX)) 
      ENDDO 
c cross-section for cold electron gas
      DO 20 IX=1,II
      X=XX(IX)
c cross-section for cold electron gas
      CALL QSN(X,0,CS,EPS) 
c total scattering cross-section
      CS=CS*(1D0+0.1D0*2D0) 
      ABSCON=dble(sigmaneu(ix,xr,ii,ab_,ab0,xnor))
      DLA=CS/(CS+ABSCON) 
      DLAMB(IX)=DLA
   20 CONTINUE
      RETURN
      END
c**********************************************************************
      SUBROUTINE REFLMATR
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NRG=7,NFEAB=3,NCA=3,ND=1)
      PARAMETER(NRAY=12,KKR=NCA*NRAY,LRAY=ND*KKR,LD=ND*NCA)
      COMMON/REFEAB/FEABREF(NFEAB)
      COMMON/RAYMATR/RGRAY(LRAY,LD),DLAMRAY(NRAY),EPSRAY(NRAY)
      DIMENSION DSTE(NRG),DS2(NRG) 
      DATA DSTE/0.5,1.,2.,4.,8.,16.,32./
      DATA DS2/0.5,0.125,0.25,0.5,1.,2.,4./
c      DATA DSTE/0.5,1.,2.,4./ 
c      DATA DS2/0.5,0.125,0.25,0.5/ 

      character*72 absfil
      character*40 filename
      character*13 info
      integer lenact,ilun,ios
      include '../spectral/inc/sitedef.inc'
c initial values to reflection matrices
      DO 200 L=1,LRAY
      DO 200 L1=1,LD
      RGRAY(L,L1)=0D0
  200 CONTINUE

      call getlun(ilun)
      absfil= XSPECDR(:lenact(XSPECDR))//'rel3.dat'
      call openwr(ilun, absfil, 'old', ' ', ' ', 0, 0, ios)
c      ilun=21
c      filename='/ymir_1/juri/PROG/REF.DAT/rel3.dat'
c      OPEN(ilun,file=filename,form='formatted',status='old')

      READ(ilun,3331) NRAY1,NC1
      DO 420 N=1,ND
      DO 420 N1=1,ND
      DO 420 I=1,NRAY
      I1=I
      READ(ilun,3332) DLA
      DLAMRAY(I)=DLA
      EPSRAY(I)=DLA/(1D0+DSQRT(1D0-DLA))**2
      DO 420 J=1,NCA
      K= J+(I-1)*NCA
      L= K+(N-1)*KKR
       READ(ilun,3333) (RGRAY(L,J1+(N1-1)*NCA),J1=1,NCA)
  420 CONTINUE
      close(ilun)
      call frelun(ilun)

 3331 FORMAT(1x,2I5)
 3332 FORMAT(1x,D13.6)
 3333 FORMAT(1x,7D13.6)
C COMPTON PART OF THE REFLECTION MATRIX
c      ilun=27
c for five Fe abundances 0.1; 0.3; 1; 3; 10.
      do ife=1,nfeab
      if(ife.eq.1) filename='ion_fe01.dat'
c      if(ife.eq.2) filename='ion_fe03.dat'
      if(ife.eq.2) filename='ion_fe1.dat'
c      if(ife.eq.4) filename='ion_fe3.dat'
      if(ife.eq.3) filename='ion_fe10.dat' 
c      if(ife.eq.1) filename='/ymir_1/juri/PROG/REF.DAT/ion_fe01.dat'
c      if(ife.eq.2) filename='/ymir_1/juri/PROG/REF.DAT/ion_fe03.dat'
c      if(ife.eq.2) filename='/ymir_1/juri/PROG/REF.DAT/ion_fe1.dat'
c      if(ife.eq.4) filename='/ymir_1/juri/PROG/REF.DAT/ion_fe3.dat'
c      if(ife.eq.3) filename='/ymir_1/juri/PROG/REF.DAT/ion_fe10.dat'
      call getlun(ilun)
      absfil= XSPECDR(:lenact(XSPECDR))//filename
      call openwr(ilun, absfil, 'old', ' ', ' ', 0, 0, ios)
c      OPEN(ilun,file=filename,form='formatted',status='old')

      READ(ilun,3330) INFO,FEABREF(ife)
      FEABREF(ife)=DLOG10(FEABREF(ife))
      do ids=1,nrg 
      CALL FREQINT(DSTE(ids),DS2(ids),ids,ilun,ife)
      enddo 
      close(ilun)
      call frelun(ilun)
c loop in Fe abundance
      enddo
 3330 FORMAT(A13,F4.1)

      RETURN
      END
c**************************************************

      SUBROUTINE FREQINT(DSTE,DSTE2,IRG,ilun,ife)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NRG=7,NFEAB=3,NCA=3,ND=1)
      PARAMETER(NCSM=50,KKC=NCA*NCSM,LCSM=ND*KKC)
      PARAMETER(IPE=NCSM+NCSM*(NCSM-1)/2,KPE=NCA*IPE)
      COMMON/RCSMA/RGCSM(KPE,NCA,NRG,NFEAB),XCSM(NCSM,NRG)
      COMMON/CIN/UANIN(NCA),AANIN(NCA)
      
c initial values to reflection matrix
      DO 201 L=1,KPE
      DO 201 L1=1,NCA
      RGCSM(L,L1,IRG,IFE)=0D0
  201 CONTINUE
c**************************************************
      READ(ilun,3331) NCSM12,NC12
      DO 421 I=1,NCSM
       READ(ilun,3332) XCSM(I,IRG)
      DO 421 I1=I,NCSM
      IP=I+(I1-1)*I1/2
      INKK=(IP-1)*NCA
      DO 421 J1=1,NCA
      L1= J1
       READ(ilun,3333) (RGCSM(J+INKK,L1,IRG,IFE),J=1,NCA)
  421 CONTINUE

 3331 FORMAT(1x,2I5)
 3332 FORMAT(1x,D13.6)
 3333 FORMAT(1x,7D13.6)
c multiplication with weights and correction for intensity
      DO 521 I=1,NCSM
      DO 521 I1=I,NCSM
      IP=I+(I1-1)*I1/2
      DFNJ=DFLOAT(NCSM-I1)+DSTE2
      XWEI=DSTE/(DFNJ*DFNJ)
      DO 521 J=1,NCA
      K=J+(IP-1)*NCA
      DO 521 J1=1,NCA
      RGCSM(K,J1,IRG,IFE)=RGCSM(K,J1,IRG,IFE)*XWEI*AANIN(J1)*
     * XCSM(I,IRG)/XCSM(I1,IRG)
  521 CONTINUE

      RETURN
      END
c**************************************************

c**************************************************
      SUBROUTINE REFLCOMP(FIN,FOUT,IFE,IRG)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NRG=7,NFEAB=3,NCA=3,ND=1)
      PARAMETER(NCSM=50,KKC=NCA*NCSM,LCSM=ND*KKC)
      PARAMETER(IPE=NCSM+NCSM*(NCSM-1)/2,KPE=NCA*IPE)

      COMMON/RCSMA/RGCSM(KPE,NCA,NRG,NFEAB),XCSM(NCSM,NRG)
      COMMON/REFEAB/FEABREF(NFEAB)
      DIMENSION FIN(KKC),FOUT(KKC),FINT(NFEAB)
      DATA CLGMIN/-30d0/,CMIN/1d-30/

      DO 30 I=1,NCSM
      DO 30 J=1,NCA
      K=J+(I-1)*NCA
c      do ire=1,nfeab
      do ire=IFE,IFE+1
      SUM=0D0
      DO 20 I1=I,NCSM
      IP=I+(I1-1)*I1/2
      KP=J+(IP-1)*NCA
      SUI=0D0
      DO 10 J1=1,NCA
      K1=J1+(I1-1)*NCA
c      SUI=SUI+RGCSM(KP,J1,IRG,ire)*FIN(K1)
      IF(ire.gt.0.and.ire.le.NFEAB)
     & SUI=SUI+RGCSM(KP,J1,IRG,ire)*FIN(K1)
   10 CONTINUE
      SUM=SUM+SUI
   20 CONTINUE
      IF(ire.gt.0.and.ire.le.NFEAB) THEN
      IF(SUM.lt.CMIN) THEN
      FINT(ire)=CLGMIN
      ELSE
      FINT(ire)=DLOG10(SUM)
      ENDIF
      ENDIF
      enddo

c interpolation to a given Fe abundance
      IF(IFE.eq.0) FOUT(K)=10D0**FINT(1)
      IF(IFE.eq.NFEAB) FOUT(K)=10D0**FINT(NFEAB)
      IF(IFE.gt.0.and.IFE.lt.NFEAB) THEN
      CALL DLINEAR(FEABREF(IFE),FINT(IFE),FEABREF(IFE+1),FINT(IFE+1),
     *  FEABLG,FLY)
c      CALL POLINTQ(FEABREF,FINT,NFEAB,FEABLG,FLY,DELY)
c      CALL RATINTQ(FEABREF,FINT,NFEAB,FEABLG,FLY,DELY)
      FOUT(K)=10D0**FLY
      ENDIF
   30 CONTINUE

      RETURN
      END
c**********************************************************************
      SUBROUTINE RAYREFL(FINT,DLA,FOUT)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NCA=3,ND=1)
      PARAMETER(NRAY=12,KKR=NCA*NRAY,LRAY=ND*KKR,LD=ND*NCA)
      COMMON/CIN/UANIN(NCA),AANIN(NCA)
      COMMON/RAYMATR/RGRAY(LRAY,LD),DLAMRAY(NRAY),EPSRAY(NRAY)
      DIMENSION  FINT(NCA),FOUT(NCA)
      DIMENSION  REFLRAY(NCA,NCA)
      save ILO

      EPSLA=DLA/(1D0+DSQRT(1D0-DLA))**2
      CALL QHUNT(DLAMRAY,NRAY,DLA,ILO)

      IF(ILO.ne.0) THEN
      EPSILO=EPSRAY(ILO)
      EPSILO1=EPSRAY(ILO+1)
      DLAILO=DLAMRAY(ILO)
c      T=(EPSLA-EPSILO)/(EPSRAY(ILO+1)-EPSILO)
      T=(DLA-DLAILO)/(DLAMRAY(ILO+1)-DLAILO)
      LAM=1
      ELSE
c if lambda < 0.01
      LAM=0
      ILO=1
      EPSILO=EPSRAY(1)
      T=0D0
      ENDIF


      DO 20 J=1,NCA
      KLO=J+(ILO-1)*NCA
      DO 20 J1=1,NCA
      IF(LAM .ne. 0) THEN
      Y1=RGRAY(KLO,J1)/EPSILO
      Y2=RGRAY(KLO+NCA,J1)/EPSILO1
      REFLRAY(J,J1)=(Y1+T*(Y2-Y1))*EPSLA
      ELSE
      REFLRAY(J,J1)=RGRAY(KLO,J1)/EPSILO*EPSLA
      ENDIF
   20 CONTINUE

      DO 30 J=1,NCA
      SUM=0D0
      DO 35 J1=1,NCA
      SUM=SUM+AANIN(J1)*FINT(J1)*REFLRAY(J,J1)
   35 CONTINUE
      FOUT(J)=SUM
   30 CONTINUE

      RETURN
      END

c**********************************************************************
