c     PROGRAM  PERMITTIVITY.F - A MEX FILE FOR MATLAB - Modified by Dijana, aug 21st, 2001
c     PRODUCES A FILE OF COMPUTED PERMITTIVITY
c     FROM THE REFLECTION COEFFICIENT FILES FOR AIR AND SAMPLE
c     ASSUMED FREQUENCY RANGE: 0.05GHz - 20.0GHz
c     ASSUMED VNA FILES SIZE: 801 POINTS
c     THE INPUT AND OUTPUT DATA ARRAYS ARE PASSED FROM/TO MATLAB
c

c     ----------------------------------------------------------------

      subroutine mexFunction(nlhs, plhs, nrhs, prhs)
      integer plhs(*), prhs(*)
      integer mxCreateFull, mxGetPr
      integer od_pr, rd_pr, pd_pr, cr_pr, a_pr

      integer nlhs, nrhs
      integer mxGetM, mxGetN, mxIsNumeric
      integer m, n, size

      real*8 open_data(801,3),refl_data(801,3)
      real*8 perm_data(801,3),cor_refl(801,3)
      real*8 a

c     check for proper number of arguments
      if (nrhs .ne. 3) then
         call mexErrMsgTxt('Three inputs required.')
      elseif (nlhs .ne. 2) then
         call mexErrMsgTxt('Two outputs required.')
      endif

c     check to see if all inputs are numeric
      if (mxIsNumeric(prhs(1)) .ne. 1) then
         call mexErrMsgTxt('Input #1 is not a numeric array')
      elseif (mxIsNumeric(prhs(2)) .ne. 1) then
         call mexErrMsgTxt('Input #2 is not a numeric array')
      elseif (mxIsNumeric(prhs(3)) .ne. 1) then
         call mexErrMsgTxt('Input #3 is not a numeric')
      endif

c     check that input #3 is scalar
      m = mxGetM(prhs(3))
      n = mxGetN(prhs(3))
      if (n .ne. 1 .or. m .ne. 1) then
         call mexErrMsgTxt('Input #3 is not a scalar')
      endif

c     get the size of the input matrices
      m = mxGetM(prhs(1))
      n = mxGetN(prhs(1))
      size = m*n

c     create matrices for the return arguments
      plhs(1) = mxCreateFull(m,n,mxREAL)
      plhs(2) = mxCreateFull(m,n,mxREAL)
      
      od_pr = mxGetPr(prhs(1))
      rd_pr = mxGetPr(prhs(2))
      a_pr = mxGetPr(prhs(3))
      pd_pr = mxGetPr(plhs(1))
      cr_pr = mxGetPr(plhs(2))

c     load the data into Fortran arrays
      call mxCopyPtrToReal8(od_pr, open_data, size)
      call mxCopyPtrToReal8(rd_pr, refl_data, size)
      call mxCopyPtrToReal8(a_pr, a, 1)
      call mxCopyPtrToReal8(pd_pr, perm_data, size)
      call mxCopyPtrToReal8(cr_pr, cor_refl, size)

c
c     call the computational routine
c
      call permCalc(open_data,refl_data,a,perm_data,cor_refl,m,n)
c
c     load the output into MATLAB arrays
c
      call mxCopyReal8toPtr(perm_data, pd_pr, size)
      call mxCopyReal8toPtr(cor_refl, cr_pr, size)

      return
      end


c     -------------------------------------------------------------------

      subroutine permCalc(open_data,refl_data,ar,perm_data,cor_refl,m,n)
      implicit none
      include 'permit5.inc'
      real*8 open_data(lotpts,3),refl_data(lotpts,3)
      real*8 perm_data(lotpts,3),cor_refl(lotpts,3)
      real*8 ar
      integer m,n
      integer i
      real*8 outlo,outhi
      real*8 openmag(fewpts),openpha(fewpts)
      real*8 sampmag(fewpts),samppha(fewpts)
      real*8 airmag(fewpts),airpha(fewpts)
      real*8 outmag(fewpts),outpha(fewpts)
c
c
c     CREATE ARRAY OF EQUALLY SPACED FREQUENCY POINTS AT WHICH 
c     MEASUREMENT DATA IS INTERPOLATED
c
c
      numpts = fewpts
      outlo = 0.05
      outhi = 20.0     
      do i = 1, numpts 
         frequency(i) = outlo + (outhi-outlo) * 
     &      (float(i-1)/float(numpts-1))
      enddo

c
c     INTERPOLATE 41 POINTS FROM 48 VNA POINTS FOR INPUT ARRAYS
c
      call interpol(open_data,openmag,openpha)
c     
c     
      call interpol(refl_data,sampmag,samppha)

c     
c     CALCULATE THEORETICAL REFLECTION COEFF FOR THE GIVEN CONDUCTOR RADIUS
c
      call refl(airmag,airpha,ar)

c
c     DO NORMALIZATION 
c
      do i = 1, numpts
         outmag(i) = sampmag(i) * airmag(i) / openmag(i)
         outpha(i) = samppha(i) + airpha(i) - openpha(i)
         if (outpha(i) .ge. 180.0) outpha(i) = outpha(i) - 360.0
         if (outpha(i) .le. -180.0) outpha(i) = outpha(i) + 360.0
      enddo
c
c     added by M.Okoniewski ( May 16, 94)
c     print out corrected values of refl. coef.
c     -------------------------------------------------begin
      do i=1, numpts
         cor_refl(i,1) = frequency(i)
         cor_refl(i,2) = outmag(i)
         cor_refl(i,3) = outpha(i)
      enddo
c
c    -----------------------------------------------------end
c
      call permsub(outmag,outpha,ar)
c
c     SAVE RESULT IN OUTPUT MATRIX
c
      do i = 1, numpts
         perm_data(i,1) = frequency(i)
         perm_data(i,2) = outmag(i)
         perm_data(i,3) = outpha(i)
      enddo
c
      return
      end
c
c     SUBROUTINE INTERPOL - GIVEN A LIST OF ARGUMENTS AND FUNCTION VALUE 
c     PAIRS, PRODUCES A SECOND LIST WITH GIVEN UPPER AND LOWER ARGUMENTS
c     AND NUMBER OF POINTS.  THE ELEMENTS IN THE SECOND LIST ARE FOUND
c     USING LINEAR INTERPOLATION BETWEEN THE APPROPRIATE PAIRS IN THE
c     INPUT LIST
c
c      
      subroutine interpol(infile,outfn1,outfn2)
      implicit none
      include 'permit5.inc'
      integer i,j,x1,x2
      real*8 inarg(lotpts), infn1(lotpts), infn2(lotpts)
      real*8 outfn1(fewpts), outfn2(fewpts)
      real*8 c
      real*8 infile(lotpts,3)
c
c     READ FREQUENCY, MAGNITUDE AND PHASE OF REFL. COEFF. FROM DATA FILE
c     INTO ARRAYS inarg(), infn1(), infn2(); 
c
c
      do i = 1, lotpts
         inarg(i) = infile(i,1)
         infn1(i) = infile(i,2)
         infn2(i) = infile(i,3)
      enddo
c
c     INTERPOLATE REFLECTION COEFFICIENT AT POINTS IN frequency() ARRAY
c
      do i = 1, numpts 
         do j = 2, lotpts
            if ( (inarg(j-1) .le. frequency(i)) .and.
     &           (inarg(j) .ge. frequency(i)) ) then 
               x1 = j-1
               x2 = j
            endif
         enddo
         c = (frequency(i) - inarg(x1)) / (inarg(x2) - inarg(x1))
c         
c        Account for a +/- 180 phase shift in phase while
c        interpolating.
c
         if ((infn2(x1) .lt. 0.0) .and. (infn2(x2) .gt. 0.0))         
     &      infn2(x1) = infn2(x1) + 360.0
         if ((infn2(x1) .gt. 0.0) .and. (infn2(x2) .lt. 0.0))
     &      infn2(x2) = infn2(x2) + 360.0
c         
c        Interpolate
c
         outfn1(i) = infn1(x1) * (1-c) + infn1(x2) * c
         outfn2(i) = infn2(x1) * (1-c) + infn2(x2) * c
      enddo
      return
      end
c      
c     REFL - GIVEN A FILE OF FREQUENCIES AND PERMITTIVITIES, THIS
c     FUNCTION APPLIES THE RATIONAL FUNCTION ADMITTANCE MODEL AND PRODUCES
c     A FILE OF REFLECTION COEFFICIENTS THAT WOULD BE MEASURED WITH AN
c     IDEAL NETWORK ANALYZER.     
c
      subroutine refl(refmag,refphase,a)     
      implicit none
      include 'permit5.inc'
      real*8 a
      integer i, j, k
      real*8 pi
      real*8 alpha(6,8),beta(6,0:8)
      real*8 refmag(fewpts),refphase(fewpts)
      real*8 mag, phase
      complex*16 nfreq(fewpts)
      complex*16 bp(0:8),cq(0:8)
      complex*16 polytop,polybot,nadmit
      complex*16 permit(fewpts),reflec(fewpts)
      data alpha/
     & -3.3749355d-01,-4.0103217d-01,-1.3372753d+00,-1.7862159d+02,1.0000000d+02,1.0000000d+03,
     & 3.4860301d+00,2.5928213d+00,-1.3626561d+02,-1.2385086d+02,1.0000000d+04,1.0000000d+04,
     & -8.7102208d-01,2.5165180d+01,2.8972912d+02,4.6871356d+02,-1.0000000d+03,1.0000000d+05,
     & 2.8321422d-01,-3.4359401d+00,-1.2845699d+02,-1.4987028d+02,9.9999997d+02,-1.0000000d+04,
     & -6.0030679d-02,8.1963925d-01,4.9055731d+01,1.0244217d+01,-1.0000014d+02,1.0000000d+04,
     & 8.0276887d-03,-2.7426544d-01,-6.3337514d+00,2.9789119d+01,9.9999218d+01,-9.9999999d+02,
     & -5.9822730d-04,2.1478253d-01,1.2824584d-01,-3.9638424d+00,9.9956912d+00,1.0000001d+02,
     & 1.8626232d-05,-2.9188608d-03,6.6881163d-01,2.7074702d-02,9.7475669d-01,-1.0000103d+01/
      data beta/
     & 1.2537824d+00,-6.2123942d+01,-2.4213368d+02,-1.3012969d+03,1.0000000d+04,-1.0000000d+05,
     & 9.1715183d+00,1.2909628d+02,3.6142773d+02,2.9512315d+03,1.0000000d+04,1.0000000d+06,
     & 5.3831550d-01,-4.0605550d+01,-9.5502322d+01,-3.0523926d+03,1.0000000d+04,-1.0000000d+06,
     & -3.8611863d-01,2.1255171d+01,6.8954080d+01,2.1045702d+03,-1.0000000d+04,1.0000000d+06,
     & -1.6519184d-01,-2.6105910d+00,8.1646959d+00,-8.0781496d+02,9.9999999d+02,-1.0000000d+05,
     & 2.0980311d-01,2.0968151d-01,1.6873807d+00,2.0962404d+02,-9.9999994d+02,1.0000000d+04,
     & -3.0957959d-02,2.0369036d-01,-1.5727025d-01,-2.5167471d+01,1.0000105d+02,-1.0000000d+03,
     & 3.4271407d-03,2.5485895d-02,1.3185953d+00,1.5036625d+00,-9.9919165d+00,1.0000007d+02,
     & -1.4119482d-04,-1.4383718d-03,9.7883954d-03,-1.8107495d-01,1.5933972d-01,-9.9927551e-01/
     
c     
      pi = dacos(-1.0d0)
c         
c      
c     ASSUME COMPLEX PERMITTIVITY OF AIR EQUAL 1
c
      do i = 1, numpts
         permit(i) = dcmplx(1,0)
      enddo            
c
c     GO THROUGH ALL DATA POINTS AND APPLY THE MODEL FOR
c     EACH POINT
c
      do i = 1, numpts      
c
c     FOR EACH POINT, CALCULATE THE NORMALIZED FREQUENCY FROM FREQ (HZ)
c
         nfreq(i) = dcmplx(0.0d0,(frequency(i)*2.0d0*pi*a))
c         
c     FOR EACH POINT, CALCULATE BP (EQUATION 5.8)
c
         do j = 1, 8
            bp(j) = (0.0d0,0.0d0)
            do k = 1, 6
               bp(j) = bp(j) + alpha(k,j)*(nfreq(i)**k)
            enddo
         enddo
         bp(0) = (0.0d0,0.0d0)
c
c     FOR EACH POINT, CALCULATE CQ (EQUATIONS 5.10, 5.11)
c
         do j = 0, 8
            cq(j) = (0.0d0,0.0d0)
            do k = 1, 6
               cq(j) = cq(j) + beta(k,j)*(nfreq(i)**k)
            enddo
         enddo
         cq(0) = cq(0) + (1.0d0,0.0d0)
c
c     FOR EACH POINT, PRODUCE THE ADMITTANCE (5.13)
c
         polytop = (0.0d0, 0.0d0)
         polybot = (0.0d0, 0.0d0)
         do j = 0, 8
            polytop = polytop + bp(j)*(cdsqrt(permit(i))**j)   
            polybot = polybot + cq(j)*(cdsqrt(permit(i))**j)
         enddo
         nadmit = polytop / polybot
c
c     PRODUCE THE REFLECTION COEFFICIENT FROM THE ADMITTANCE
c
         reflec(i) = (1 - nadmit)/(1 + nadmit)
      enddo
c
c
      do i = 1, numpts
      mag = cdabs(reflec(i))
      phase = datan2( dimag(reflec(i)),dreal(reflec(i)) )
      phase = phase * 57.29578
      refmag(i)=mag
      refphase(i)=phase
      enddo
      return
      end
c
c     SUBROUTINE PERMSUB
c
c
      subroutine permsub(outmag,outpha,a)
      implicit none
      include 'permit5.inc'
      real*8 a
      integer  i, j, k
      real*8 pi, phi
      real*8 alpha(6,8),beta(6,0:8)
      real*8 refmag,refpha,refreal,refimag
      real*8 outmag(fewpts),outpha(fewpts)
      complex*16 reflec
      complex*16 nfreq(fewpts)
      complex*16 bp(0:8),cq(0:8),poly(9)
      complex*16 nadmit(fewpts),permit(fewpts),roots(8),answer
      data alpha/
     & -3.3749355d-01,-4.0103217d-01,-1.3372753d+00,-1.7862159d+02,1.0000000d+02,1.0000000d+03,
     & 3.4860301d+00,2.5928213d+00,-1.3626561d+02,-1.2385086d+02,1.0000000d+04,1.0000000d+04,
     & -8.7102208d-01,2.5165180d+01,2.8972912d+02,4.6871356d+02,-1.0000000d+03,1.0000000d+05,
     & 2.8321422d-01,-3.4359401d+00,-1.2845699d+02,-1.4987028d+02,9.9999997d+02,-1.0000000d+04,
     & -6.0030679d-02,8.1963925d-01,4.9055731d+01,1.0244217d+01,-1.0000014d+02,1.0000000d+04,
     & 8.0276887d-03,-2.7426544d-01,-6.3337514d+00,2.9789119d+01,9.9999218d+01,-9.9999999d+02,
     & -5.9822730d-04,2.1478253d-01,1.2824584d-01,-3.9638424d+00,9.9956912d+00,1.0000001d+02,
     & 1.8626232d-05,-2.9188608d-03,6.6881163d-01,2.7074702d-02,9.7475669d-01,-1.0000103d+01/
      data beta/
     & 1.2537824d+00,-6.2123942d+01,-2.4213368d+02,-1.3012969d+03,1.0000000d+04,-1.0000000d+05,
     & 9.1715183d+00,1.2909628d+02,3.6142773d+02,2.9512315d+03,1.0000000d+04,1.0000000d+06,
     & 5.3831550d-01,-4.0605550d+01,-9.5502322d+01,-3.0523926d+03,1.0000000d+04,-1.0000000d+06,
     & -3.8611863d-01,2.1255171d+01,6.8954080d+01,2.1045702d+03,-1.0000000d+04,1.0000000d+06,
     & -1.6519184d-01,-2.6105910d+00,8.1646959d+00,-8.0781496d+02,9.9999999d+02,-1.0000000d+05,
     & 2.0980311d-01,2.0968151d-01,1.6873807d+00,2.0962404d+02,-9.9999994d+02,1.0000000d+04,
     & -3.0957959d-02,2.0369036d-01,-1.5727025d-01,-2.5167471d+01,1.0000105d+02,-1.0000000d+03,
     & 3.4271407d-03,2.5485895d-02,1.3185953d+00,1.5036625d+00,-9.9919165d+00,1.0000007d+02,
     & -1.4119482d-04,-1.4383718d-03,9.7883954d-03,-1.8107495d-01,1.5933972d-01,-9.9927551e-01/
    
c     
      pi = dacos(-1.0d0)
c
c     READ MAGNITUDE AND PHASE OF REFLECTION COEFFICIENT,
c     AND CONVERT TO NORMALIZED ADMITTANCE
c
      do i = 1, numpts
         refmag=outmag(i)
         refpha=outpha(i) 
         refpha = refpha/57.2957795131
         refreal = refmag * dcos(refpha)
         refimag = refmag * dsin(refpha)
         reflec = dcmplx(refreal,refimag)
         nadmit(i) = (1 - reflec)/(1 + reflec)
      enddo
c
c     GO THROUGH ALL DATA POINTS AND CALCULATE THE PERMITTIVITY FOR
c     EACH POINT
c
      do i = 1, numpts      
c
c     FOR EACH POINT, CALCULATE THE NORMALIZED FREQUENCY FROM FREQ (HZ)
c
         nfreq(i) = dcmplx(0.0d0,(frequency(i)*2.0d0*pi*a))
c         
c     FOR EACH POINT, CALCULATE BP (EQUATION 5.8)
c
         do j = 1, 8
            bp(j) = (0.0d0,0.0d0)
            do k = 1, 6
               bp(j) = bp(j) + alpha(k,j)*(nfreq(i)**k)
            enddo
         enddo
         bp(0) = (0.0d0,0.0d0)
c
c     FOR EACH POINT, CALCULATE CQ (EQUATIONS 5.10, 5.11)
c
         do j = 0, 8
            cq(j) = (0.0d0,0.0d0)
            do k = 1, 6
               cq(j) = cq(j) + beta(k,j)*(nfreq(i)**k)
            enddo
         enddo
         cq(0) = cq(0) + (1.0d0,0.0d0)
c
c     FOR EACH POINT, PRODUCE THE COEFFICIENTS OF THE POLYNOMIAL (5.14)
c
         do j = 1, 9
            poly(j) = bp(j-1) - nadmit(i)*cq(j-1)            
         enddo
c
c     FIND THE PHYSICALLY MEANINGFUL ROOTS OF THE POLYNOMIAL
c
         call zroots(poly,8,roots,.true.)
         do j = 1, 8
            phi = datan2(dimag(roots(j)),dreal(roots(j)))
            phi = phi * 57.2958
            if(phi .lt. 1.0d-1 .and. phi .gt. -45.0d0 .and.
     +         cdabs(roots(j)) .le. 10.0) then
                  answer = roots(j)
            endif
         enddo
c
c     WRITE ANSWER TO PERMITTIVITY ARRAY
c
         permit(i) = answer*answer
      enddo
c
      do i = 1, numpts
         outmag(i)= dreal(permit(i))
         outpha(i)= -dimag(permit(i)) 
      enddo      
      return
      end
c      
c
      subroutine zroots(a,m,roots,polish)
c
c     given the degree M and the M+1 complex coeeficients A of the polynomial
c
c               A(1)+A(2)x+A(3)x**2+.....+A(M+1)x**(M)
c
c     this routine successively calls LAGUER and finds all M complex ROOTS. 
c     The logical variable POLISH should be input as .TRUE. if polishing
c     by Laguerre's method is desired, .FALSE. if polishing will subsequently 
c     be done by other means.
c
      implicit none
      real*8 eps,maxm
      parameter( eps=3.0d-17,maxm=101 )
      integer*4 m
      complex*16 a(m+1),roots(m),ad(maxm),x,b,c
      integer*4 j,jj,i
      logical polish
c
      do j=1,m+1
         ad(j)=a(j)
      enddo
      do j=m,1,-1
         x=(0.0d0,0.0d0)
         call laguer(ad,j,x,eps,.false.)
         if(dabs(dimag(x)) .le. 2.0d0*eps**2*dabs(dreal(x))) 
     &      x=dcmplx(dreal(x),0.0d0)
         roots(j)=x
         b=ad(j+1)
         do jj=j,1,-1
            c=ad(jj)
            ad(jj)=b
            b=x*b+c
         enddo
      enddo
      if(polish) then
         do j=1,m
            call laguer(a,m,roots(j),eps,.true.)
         enddo
      endif
      do j=2,m
         x=roots(j)
         do i=j-1,1,-1
         if(dreal(roots(i)) .le. dreal(x)) goto 10
         roots(i+1)=roots(i)
      enddo
      i=0
10    roots(i+1)=x
      enddo
      return
      end
c
c
      subroutine laguer(a,m,x,eps,polish)
c
c     Given the degree M and the M+1 complex coefficients A of the polynomial
c
c         A(1)+A(2)x+a(3)x**2+....+A(M+1)x**M
c
c     and given the desired fractional accuracy EPS, and given a complex value
c     X, this routine improves X by Laguerre's method until it converges to a 
c     root of the given polynomial.  For normal use POLISH should be input as
c    .FALSE..  When POLISH is input as .TRUE., the routine ignores EPS and 
c     instead attempts to improve X (assumed to be a good initial guess to the
c     desired root) to the achievable roundoff limit.
c
      implicit none
      real*8 epss,maxit
      parameter( epss=3.0d-17,maxit=100 )
      integer*4 m
      complex*16 A(m+1),x,dx,x1,b,d,f,g,h,sq,gp,gm,g2,zero
      integer*4 iter,j
      real*8 dxold,err,abx,cdx,eps
      logical polish
c
      zero=(0.0d0,0.0d0)
      dxold=cdabs(x)
      do iter=1,maxit
         b=a(m+1)
         err=cdabs(b)
         d=zero
         f=zero
         abx=cdabs(x)
         do j=m,1,-1
            f=x*f+d
            d=x*d+b
            b=x*b+a(j)
            err=cdabs(b)+abx*err
         enddo
         err=epss*err
         if(cdabs(b) .le. err) then
            dx=zero
            return
         else
            g=d/b
            g2=g*g
            h=g2-2.0d0*f/b
            sq=dsqrt(dfloat(m-1))*cdsqrt((m*h-g2))
            gp=g+sq
            gm=g-sq
            if(cdabs(gp) .lt. cdabs(gm)) gp=gm
            dx=m/gp
         endif
         x1=x-dx
         if(x .eq. x1) return
         x=x1
         cdx=cdabs(dx)
         if(iter .gt. 6 .and. cdx .ge. dxold) return
         dxold=cdx
         if( .not. polish) then
             if(cdx .le. eps*cdabs(x)) return
         endif
      enddo
      pause ' too many iterations'
      return
      end
