program trinve
******************************************************************************
*     this program inverts transmission data of a thin film on a plan-parallel
*     substrate. On input it wants to know precisely the substrate optical
*     constant, the thickness, and the thickness of the film.
*     The program samples segments of one period of Fabry-Perod oscillations
*     and fits the thin-film-on-plan-parallel-substrate formula to this.
******************************************************************************
      parameter(num=100000)
      real x(num),y(num),ylo,yhi
      real xsb(num),dsub,pl,pr,psil,psir,xl,xr,p,psi2,pi
      real xav,ymn,ymx,dfilm,kd
      integer i,i1,i2,imx,ihi(num),ilo(num),nmx,nperiod,jlo,jhi,nsub
      complex nsb(num),znw
      character*40 flin,flout,flisub
      pi=4.*atan(1.)
      write(*,*) 'give inputfilename '
      read(*,'(a40)') flin
      write(*,*) 'give inputfilename containing 
     * substrate dielectric function'
      read(*,'(a40)') flisub
      write(*,*) 'give substrate thickness '
      read(*,*) dsub
      write(*,*) 'give film thickness '
      read(*,*) dfilm
      open(14,file=flisub)
      do 5 i=1,num
       read(14,*,END=6) xsb(i),e1,e2
       nsb(i)=csqrt(cmplx(e1,e2))
5     continue
6     nsub=i-1
      close(14)       
      open(14,file=flin)
      write(*,*) 'give outputfilename '
      read(*,'(a40)') flout
      open(23,file=flout)
      do 100 i=1,num
       read(14,*,END=101) x(i),y(i)
       y(i)=1/y(i)
100   continue
101   imx=i-1 
      close(14)
      write(*,*) 'how many measurement points per period ? '
      read(*,*) nperiod
      n=1
      i1=1
150   ylo=1e8
      yhi=0
      do 200 i=i1,i1+0.75*nperiod
*      find nth maximum
c       write(*,*) yhi,ylo,y(i)
       if (y(i).gt.yhi) then
        ihi(n)=i
        yhi=y(i)
       endif
200   continue
      i1=ihi(n)
      do 205 i=i1,i1+1.25*nperiod
*      find nth minimum
       if (y(i).lt.ylo) then
        ilo(n)=i
        ylo=y(i)
       endif
205   continue
      i1=ilo(n)
      if (.not.(n.eq.1)) then
       if (n.gt.5) then
        nperiod=(ihi(n)-ihi(n-5))/5
       else
        nperiod=(ihi(n)-ihi(1))/(n-1)
       endif
      endif  
      write(*,*) 'n,ilo(n),ihi(n),period',n,ilo(n),ihi(n),nperiod
      i2=i1+1.25*nperiod
      if (i2.gt.imx) goto 151
      n=n+1
      goto 150
151   nmx=n 
      do 300 n=1,nmx-1
       jlo=ilo(n)
       jhi=ihi(n)
       ymn=y(ilo(n))
       ymx=(y(ihi(n))+y(ihi(n+1)))/2
       xav=x(ilo(n))
c      find the substrate psi2 and p for xav by interpolation
       do 250 i=1,nsub
        if (xsb(i).lt.xav) then
         il=i
         ir=i+1
        endif
250    continue
       pl=real(nsb(il))
       pr=real(nsb(ir))
       psil=dsub*2*pi*xsb(il)*aimag(nsb(il))
       psir=dsub*2*pi*xsb(ir)*aimag(nsb(ir))
       xl=xsb(il)
       xr=xsb(ir)
       p=pl+(pr-pl)*(xav-xl)/(xr-xl)        
       psi2=psil+(psir-psil)*(xav-xl)/(xr-xl)
       call filsub(p,psi2,ymx,ymn,znw)
       kd=dfilm*2*pi*xav
       write(*,*) xav,-real(znw/kd),aimag(znw/kd)
       write(23,*) xav,-real(znw/kd),aimag(znw/kd)
300   continue
      close(23)
      end
******************************************************************************
      
******************************************************************************
      subroutine filsub(p,psi2,gmx,gmn,zout)
      complex zout
      real psi2,p,chp,shp,p2,p4,gmx,gmn
     * ,gav,gam,a,b,d,at,dt,ff,denom,n,m,ynw,zz,xnw
      gav=(gmx+gmn)/2
      gam=(gmx-gmn)/2
      p2=p**2
      p4=p**4
      chp=cosh(2*psi2)
      shp=sinh(2*psi2)
      a=(1+3*p2)*chp/(4*p2)+(3+p2)*shp/(4*p)
      b=(1+p2)*chp/(8*p2)+shp/(4*p)
      d=(1+6*p2+p4)*chp/(8*p2)+(1+p2)*shp/(2*p)-gav
      at=a/b
      dt=d/b
      ff=(8*p2*gam/(p2-1))**2
c      write(*,*) 'at,dt,ff ',at,dt,ff
      denom=(at-2)**2-4*p2
      n=(ff+4*p2*dt-(p2+dt-1)**2)/denom
      m=((at-2)*(p2+dt-1)-2*p2*at)/denom
c      write(*,*) 'm,n ',m,n
      ynw=-m+sqrt(m**2+n)
      zz=-at*ynw-dt
c      write(*,*) 'at,ynw,dt,zz ',at,ynw,dt,zz
      zz=zz-ynw**2
      if (zz.lt.0) then
       xnw=0
      else 
       xnw=sqrt(zz)
      endif 
      zout=cmplx(xnw,ynw)
      return
      end
******************************************************************************