program kramerm
************************************************************************
* Kramers-Kronig analysis of reflectance data.                         *
* Designed for metals.                                                 *
* From standard input it reads:                                        *
*    inputfilename  (x y format)                                       *
*    NLOW NHIG INC                                                     *
* These are the lowest, highest and incremental integer numbers which  *
* indicate the data points for which the transform has to be made.     *
* A file containing x epsilon1 epsilon2 is flushed to standard output. *
************************************************************************
      parameter(nfys=10000)
       real x,r,lnr(0:nfys),w(0:nfys),teta(nfys),xlast,tol
       real work(0:nfys)
       complex ce,cn,cr
       integer i,ndat,nlow,nhig,inc,tel
       character*40 flin
****** Never used: phi(nfys),flsig,fleps,cea,fln,fllnr,fllos
       tol=0.000001
       read(*,'(a40)') flin
       read(*,*) nlow,nhig,inc
       open(17,file=flin)
       open(18,file='kramer.log')
       lnr(0)=0.
       w(0)=0.
       xlast=-1e20
       tel=0
       do 400 i=1,nfys
         read(17,*,END=401) x,r
         if (abs(x-xlast).lt.tol) then
          xlast=x
          goto 400
         else
          tel=tel+1
          w(tel)=x
          xlast=x
         endif
         if (r.lt.0) then
          write(18,*)  'Reflectivity ',r,' smaller than 0 at
     *   frequency ',w(tel), 'Corrective action taken.'
          lnr(tel)=lnr(tel-1)
          goto 400
         endif
         if (r.gt.1) then
          write(18,*)  'Reflectivity ',r,' larger than 1 at
     *   frequency ',w(i), 'Corrective action taken.'
          lnr(tel)=lnr(tel-1)
          goto 400
         endif
         lnr(tel)=0.5*alog(r)
400    continue
c        That was
c401    ndat=tel-1
401    ndat=tel
       if (nlow.lt.1) nlow=1
       if (nhig.ge.ndat) nhig=ndat-1
       call kkr(work,w,lnr,ndat,teta,nlow,nhig,inc)


* The phase angle \teta is known modulo \pi. This is
* due to the fact that the sign of \sqrt(R) is not fixed experimentally.
* However, the way the transform is made, guarantees that \teta=0
* for x=\infty. This implies that we must use r=(n-1)/(n+1), as this
* has a phase angle zero (not pi) for Re(n)>1, Im(n)=0. The latter
* condition is automatically satisfied in the high frequency limit
* of a physical response function.
* Furthermore we expect \teta>0 in most cases:
* Due to small numerical perturbations it may happen
* that teta obtains the wrong sign near a zero crossing. As far
* as I can see an exact transform should automatically produce the
* correct sign. The sign-correction results in
* a response function, which is no longer analytical, so it is
* awkwart. However a wrong sign of \teta is worse: It violates
* causality, as it corresponds to a negative value of the optical
* constant (and hence epsilon).
* A problem arises if we work with pseudo dielectric functions due to
* p-polarized light: In this case negative values of \phi ARE
* physical. I don't know whether a cure exists for all cases. The best
* thing is to check the resulting phase function for possible
* regions of negative phi. Hence:
       do 600 i=nlow,nhig,inc
* depending on the kind of reflection spectrum one may
* decide to do this:  if (teta(i).lt.0) teta(i)=-teta(i)
        x=w(i)
        cr=exp(cmplx(lnr(i),teta(i)))
        if (cr.eq.1.) goto 600
        cn=(1.+cr)/(1.-cr)
        ce=cn**2
        write(*,*) x,real(ce),aimag(ce)
600    continue
       close(18)
      END
***********************************************************************


************************************************************************
***** The input file y is assumed to be symmetric in frequency.        *
***** The transformed output file z is then asymmetric in frequency.   *
***** z(j)=-\frac{2j}{\pi}\int_0^{\infty}\frac{y(i)-y(j)}{i^2-j^2}di   *
************************************************************************
      SUBROUTINE kkr(wrk,x,y,ntot,z,nlow,nhig,inc)
      REAL x(0:ntot),y(0:ntot),z(ntot),wrk(0:ntot),pi,sum,klad
      INTEGER i,j
      pi=4.*atan(1.)
      DO 10 j=nlow,nhig,inc
       DO 20 i=0,ntot
        IF (i.NE.j) wrk(i)=x(j)*(y(i)-y(j))/(x(i)*x(i)-x(j)*x(j))
20     CONTINUE
        wrk(j)=(wrk(j-1)+wrk(j+1))/2.
       sum=0.
       DO 30 i=1,ntot-1
        sum=sum+((wrk(i)+wrk(i+1))/2.)*(x(i+1)-x(i))
30     continue
****** There are two ways to calculate contribution from large $\omega$.
****** 1:
****** We assume that $y$ and $z$ are constant for all $\omega >$ ntot.
****** This happens for large $\omega$, where $\eps=eps_\infty.
****** We use an analytical extension to calculate the contribution
****** to z in this frequency range.
       klad=(x(ntot)-x(j))/(x(ntot)+x(j))
       if (klad.le.0) write(18,*) 'warning ',j,x(j),x(ntot),klad
       sum=sum-0.5*(y(ntot)-y(j))*alog(klad)
       z(j)=-(2./pi)*sum
******
****** 2:
****** We assume that $y$ is varying inversely to the forth degree 
****** of $\omega$ for $\omega >$ ntot.
****** In this case we should add another term to the previous integral.
****** See F.Wooten, Optical properties of solids, p.249.
c       do 40 i=0,10
c        z(j)=z(j)+(4./pi)*((x(j)/x(ntot))**(2*i+1))/((2.*i+1.)**2)
c40     continue
10    CONTINUE
      RETURN
      END
***********************************************************************