C
C   THIS PROGRAM CALIBRATES THE SCATTEROMETERS FOR
C   UNIFORMLY DISTRIBUTED TARGET.             10-17-91 
C   REVISED FOR NCRB SPHERE DATA              12-31-91
C   REVISED FOR CHAMBER SPHERE DATA            4-25-92
C   CORRECT REFERENCE DATA W.R.T. CHAMBER DATA   5-4-92
C                                                       

      subroutine subcalc(nntrace,ttype,aang,llabel,ttitle)
	

      complex s0(21),esph0(4,21),esphd(15,15,4,21)
      complex esoil(100,4,21)
      complex scal(100,4,21),d(15,15,21,4,4)
      complex b(16,16,21),y(16,21),x(16,21),mueller(4,4)
      complex z1,z2,cj,muel(21,4,4),yy(16,21)
      real sig(4,21),sigav(4),alpha,eta,rsphfld,rsphcham
      real angle,fstart,fstop,pi,range,height,correl
      real zetmin,zetmax,psimin,psimax
      real alp(21),zet(21),alpx(21),zetx(21)
      integer ntrace,ndata,np,nzeta,npsi,nskip,nntrace
      character*30 reffile,scalfile,smeafile,skyfile,gndfile
      character*30 msoutfile,fasoutfile,fasgrfile
      character aang*2,ttype*1,llabel*11,ttitle*4
      character ang*2,type*1,label*11,title*4

      common /sphere/ s0,esph0,esphd
      common /target/ esoil
      common /distor/ d
      common /calmat/ b  
      common /measur/ y,yy
      common /xcorr/ x   
      common /muell/ mueller,muel
      common /scatt/ sig,sigav      
      common /phaze/ alpha,eta,alphax,etax
      common /fas/ alp,zet,alpx,zetx
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /para2/ range,height,rsphfld,rsphcham
      common /para3/ nzeta,npsi,zetmin,zetmax,psimin,psimax
      common /file1/ msoutfile,fasoutfile,gndfile,reffile,smeafile,it
      common /names/ ang,type,label,title

                                
        ntrace=nntrace
	type=ttype
	ang=aang
	label=llabel
	title=ttitle

C  READ INPUT DATA
C                 
      call readinc
c
c                            print *,'read data'
C
C  CALCULATE MATRIX "D" (calibration matrix)
C
c                            print *,'going to dmatrixc'                          
      call dmatrixc                             
c                            print *,'calculated D matrix'
C
C  CALCULATE MATRIX "B" (cross-correlation calibration matrix)
C
      call bmatrixc                                       
c                            print *,'calculated B matrix'
       do 5 i=1,16
       do 5 id=1,ndata
          x(i,id) = (0.0,0.0)
          yy(i,id) = (0.0,0.0)
   5  continue
                          
      do 10 it = 1,ntrace
C        print *,'it=',it
C
C  CALCULATE MEASURED VECTOR "Y" (measured cross-correlation vector)
C                     
      call yvectorc(it)
C
C  SOLVE LINEAR EQUATION (Y = B X, where X is unknown.)
C                       
      call solveqc 
c   10 continue    
c                            print *,'calculated X matrix'
C
C  FIND BACKSCATTERING COEFFICIENTS THRU MUELLER MATRIX
C
      call sigmac   
c                            print *,'calculated Mueller matrix'
C
C  PRINT OUT
C
      call prtoutc   
C
10    continue
      return
      end
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine readinc
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C        
C  FIELD DATA POL. SEQ. = VH - VV - HH - HV      
C  CHAMB DATA POL. SEQ. = VH - VV - HH - HV      
C
C  STORED POL. SEQ.     = VV - VH - HV - HH
C                                            
      complex esky(4,21)
      complex s0(21),esph0(4,21),esphd(15,15,4,21)
      complex esoil(100,4,21)
      complex cj,temp                                       
      real angle,fstart,fstop,pi,range,height,correl
      real zetmin,zetmax,psimin,psimax,rsphfld,rsphcham
      integer ntrace,ndata,np,nzeta,npsi
      character*30 reffile,scalfile,smeafile,skyfile,gndfile
      character*30 msoutfile,fasoutfile,fasgrfile
      character ang*2,type*1,label*11,title*4

      common /sphere/ s0,esph0,esphd
      common /target/ esoil
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /para2/ range,height,rsphfld,rsphcham
      common /para3/ nzeta,npsi,zetmin,zetmax,psimin,psimax
      common /file1/ msoutfile,fasoutfile,gndfile,reffile,smeafile,it
      common /names/ ang,type,label,title
                         
      scalfile='../util.lib/sph_14_c'
      reffile=title//'_ref.'//type
      smeafile='../util.lib/sph_c_cham'
      skyfile =title//'_'//ang//'s.'//type
      gndfile =type//label       
      msoutfile='ac'//label 
      fasoutfile='pc'//label
      
      if (ang .eq. '20') then
      angle = 20.0
      else if(ang .eq. '30') then
      angle = 30.0
      else if(ang .eq. '40') then
      angle = 40.0
      else if(ang .eq. '50') then
      angle = 50.0
      else if(ang .eq. '60') then
      angle = 60.0
      else
      angle = 70.0
      endif


      fstart=5.50
      fstop=5.10
      ntrace=ntrace
      npts=21
      ndata=21
      nskip=2
      range=12.
      offset=1.15
      rsphfld=10.85
      rsphcham=9.15
      nzeta=15.
      zetmin=-10.5
      zetmax=10.5

C                 
C  INITIAL PARAMETER
C
      pi = 4.*atan(1.)
      cj = cmplx(0.0,1.0)   
      np = 4            
      range = range - offset
      height = range * cos(angle*pi/180.)

      npsi = nzeta 
      zetmax =zetmax *pi/180.
      zetmin =zetmin *pi/180.
      psimin = zetmin 
      psimax = zetmax         
      print *,ntrace, ndata, height,range,rsphcham
C
C  OPEN FILES
C
        open(10,file=scalfile)
        open(11,file=reffile)
        open(12,file=smeafile) 
        open(13,file=skyfile) 
        open(14,file=gndfile) 
C
C   READ THEORETICAL VALUE OF SPHERE
C
      do 10 i=1,npts
        read (10,*) dum,s0(i)
   10 continue     
C
C  SPHERE DATA MEASURED AT BORESIGHT ONLY (AT FIELD), AND REARRANGE THEM.
C
      do 20 ip=1,np
        read(11,100)
      do 20 id=1,npts
         read(11,*) esph0(ip,id)
   20 continue               
     
      do 25 id=1,npts
        temp        = esph0(1,id)
        esph0(1,id) = esph0(2,id)
        esph0(2,id) = temp

        temp        = esph0(3,id)
        esph0(3,id) = esph0(4,id)
        esph0(4,id) = temp
   25 continue        
C
C  SPHERE DATA MEASURED FOR ANTENNA PATTERN AT CHAMBER
C
      do 30 iz=1,nzeta
      do 30 is=1,npsi
      do 30 ip=1,np 
         read(12,100)
         read(12,*) (esphd(iz,is,ip,id),id=1,npts)
   30 continue    

      do 35 iz=1,nzeta
      do 35 is=1,npsi
      do 35 id=1,ndata
        temp              = esphd(iz,is,1,id)
        esphd(iz,is,1,id) = esphd(iz,is,2,id)
        esphd(iz,is,2,id) = temp

        temp              = esphd(iz,is,3,id)
        esphd(iz,is,3,id) = esphd(iz,is,4,id)
        esphd(iz,is,4,id) = temp
   35 continue
C
C SKY DATA   
C                  
c      if (skyfilexist) then   

       do 40 ip=1,np
         read(13,100)
       do 40 id=1,npts
          read(13,*) esky(ip,id)
   40  continue     

c      else          

c       do 45 ip=1,np
c       do 45 id=1,npts
c          esky(ip,id) = (0.0, 0.0)
c   45  continue     

c      endif
C
C TARGET (SOIL) DATA
C     
      do 50 it=1,ntrace                
      do 50 ip=1,np
        read(14,100)
      do 50 id=1,npts
         read(14,*) esoil(it,ip,id)
         esoil(it,ip,id) = esoil(it,ip,id) - esky(ip,id)
   50 continue    

      print *,'soilfile-read'    
      do 55 it= 1,ntrace
      do 55 id=1,npts
        temp           = esoil(it,1,id)
        esoil(it,1,id) = esoil(it,2,id)
        esoil(it,2,id) = temp

        temp           = esoil(it,3,id)
        esoil(it,3,id) = esoil(it,4,id)
        esoil(it,4,id) = temp
   55 continue        
  100 format(1x) 

      do i=1,5 
        close(9+i)
      enddo  
C
C  OPEN FILES
C
        open(20,file=msoutfile)
        open(21,file=fasoutfile)
C
      return                                                              
      end
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine dmatrixc
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                                                  
      complex s0(21),esph0(4,21),esphd(15,15,4,21),dvv,dhh
      complex esoil(100,4,21),cj,cch,cfi(21)
      complex d(15,15,21,4,4),a,c,alph,beta,rt,cdat(15,15,21)
      complex frt,falph,fbeta,facrt(21),facalp(21),facbet(21)
      complex c1,c2,c3,c4,ac,ach,afi

      common /sphere/ s0,esph0,esphd
      common /distor/ d
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /para2/ range,height,rsphfld,rsphcham
      common /para3/ nzeta,npsi,zetmin,zetmax,psimin,psimax
       
C
      nzb = (nzeta+1)/2
      npb = (npsi+1)/2                
C
C  CORRECT REFERENCE DATA W.R.T. CHAMBER DATA
C
c      ndt = ndata-2*nskip
c
c      fhhvv = 0.0
c      fvhvv = 0.0
c      fhvvv = 0.0
c      chhvv = 0.0
c      cvhvv = 0.0
c      chvvv = 0.0
c      do 2 id=1+nskip,ndata-nskip
c
c       fhhvv = fhhvv + cabs(esph0(4,id)/esph0(1,id))/ndt
c       fvhvv = fvhvv + cabs(esph0(2,id)/esph0(1,id))/ndt
c       fhvvv = fhvvv + cabs(esph0(3,id)/esph0(1,id))/ndt
c       chhvv =chhvv+cabs(esphd(nzb,npb,4,id)/esphd(nzb,npb,1,id))/ndt
c       cvhvv =cvhvv+cabs(esphd(nzb,npb,2,id)/esphd(nzb,npb,1,id))/ndt
c       chvvv =chvvv+cabs(esphd(nzb,npb,3,id)/esphd(nzb,npb,1,id))/ndt
c   2  continue
c        print *,'fhhvv,fvhvv,c...=',fhhvv,fvhvv,chhvv,cvhvv
c        print *,'esph0(4,11)=',20.*alog10(cabs(esph0(4,11)))
c      do 4 id=1,ndata
c        esph0(4,id) = esph0(4,id) *chhvv/fhhvv
c        esph0(2,id) = esph0(2,id) *cvhvv/fvhvv
c        esph0(3,id) = esph0(3,id) *chvvv/fhvvv
c   4  continue
c        print *,'esph0(4,11) cor=',20.*alog10(cabs(esph0(4,11)))
C
C  FIND CORECTION FACTORS WRT  R(TILDE), T(TILDE)
C                                                 
      do 10 id=1,ndata        
 
C
C  CHAMBER MEASURED SPHERE DATA (BORE SIGHT)
C
        a = (esphd(nzb,npb,2,id)*esphd(nzb,npb,3,id))
     &         /(esphd(nzb,npb,1,id)*esphd(nzb,npb,4,id))
        c = (1.0 - csqrt(1.0-a)) / csqrt(a)
        
        print*,'chamber c =',c
        cfas=phasec(c)
        if(cfas.lt.0.0) c=-c

        print*, 'chamber c =', c, 'phase =',cfas

C
c        err =sqrt( cabs(esphd(nzb,npb,3,id))
c     &                      /cabs(esphd(nzb,npb,2,id)) )
        err=1.0
        
        rt = esphd(nzb,npb,1,id) / (1.0 + c**2) / s0(id)
        beta = 2.0*c /(1.0+c**2) *esphd(nzb,npb,4,id)
     &                        /esphd(nzb,npb,2,id)   /err
        alph = (1.0+c**2) /(2.0*c) *esphd(nzb,npb,2,id)
     &                        /esphd(nzb,npb,1,id)   *err
C
C  FIELD MEASURED SPHERE DATA (BORE SIGHT)
C                                                    
        a = (esph0(2,id)*esph0(3,id)) /(esph0(1,id)*esph0(4,id))
        c = (1.0 - csqrt(1.0-a)) / csqrt(a)
C
C  CORRECT "c" VALUES
C           
        cfi(id)=c            
        cfas=phasec(c)
        if(cfas.lt.0.0) c=-c
c       print *,id,phasec(c),phasec(cfi(id))
C
        print *, 'field c =', c

        err = sqrt( cabs(esph0(3,id)) /cabs(esph0(2,id)) )
        frt = esph0(1,id) / (1.0 + c**2) / s0(id)
        fbeta = 2.0*c /(1.0+c**2) *esph0(4,id)/esph0(2,id) 
c     &                              /err
        falph = (1.0+c**2) /(2.0*c) *esph0(2,id)/esph0(1,id)
c     &                              *err
C
C  CORECTION FACTORS
C                                                           
        facrt(id) = frt / rt
        facalp(id) = falph /alph
        facbet(id) = fbeta /beta
   10 continue
C
C  CALCULATE THE ELEMENTS OF THE "D" MATRIX
C
      do 20 iz=1,nzeta
      do 20 is=1,npsi   
      do 20 id=1,ndata
        a = (esphd(iz,is,2,id)*esphd(iz,is,3,id))
     &         /(esphd(iz,is,1,id)*esphd(iz,is,4,id))
        c = (1.0 - csqrt(1.0-a)) / csqrt(a)
C
C  CORRECT "c" VALUES
C                       
        cfas=phasec(c)
        if(cfas.lt.0.0) c=-c
C
        err = sqrt( cabs(esphd(iz,is,3,id))
     &                    /cabs(esphd(iz,is,2,id)) )
        rt = esphd(iz,is,1,id) / (1.0 + c**2) / s0(id) 
        beta = 2.0*c /(1.0+c**2) *esphd(iz,is,4,id)/esphd(iz,is,2,id)
c     &                         /err
        alph = (1.0+c**2) /(2.0*c) *esphd(iz,is,2,id)/esphd(iz,is,1,id)
c     &                         *err

        rt = rt *facrt(id)                     
        beta = beta *facbet(id)
        alph = alph * facalp(id)

        d(iz,is,id,1,1) = rt * 1.0
        d(iz,is,id,1,2) = rt * c * alph
        d(iz,is,id,1,3) = rt * c
        d(iz,is,id,1,4) = rt * c**2 * alph
        d(iz,is,id,2,1) = rt * c
        d(iz,is,id,2,2) = rt * alph
        d(iz,is,id,2,3) = rt * c**2
        d(iz,is,id,2,4) = rt * c * alph
        d(iz,is,id,3,1) = rt * c * beta
        d(iz,is,id,3,2) = rt * c**2 * alph * beta
        d(iz,is,id,3,3) = rt * beta
        d(iz,is,id,3,4) = rt * c * alph * beta
        d(iz,is,id,4,1) = rt * c**2 *beta
        d(iz,is,id,4,2) = rt * c * alph * beta
        d(iz,is,id,4,3) = rt * c * beta
        d(iz,is,id,4,4) = rt * alph * beta
   20 continue      

      return
      end
C       
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
        FUNCTION PHASEC(Z)
C                                                           
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        COMPLEX Z 
        PI=4.*ATAN(1.)
        X=REAL(Z)
        Y=AIMAG(Z)
        PHASEC=(180./PI)*ATAN2(Y,X)
        RETURN
        END
C                                                                  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine bmatrixc
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                                                  
      complex b(16,16,21),sums,sumz,funb,sums1,sumz1,dsq
      complex d(15,15,21,4,4),z1,z2,cj,sumdsq(16,16,21)
      real dz,ds,dzrad,dsrad,zeta,psi,trig
      integer i,j,l,m,iz,is
      
      common /distor/ d            
      common /calmat/ b
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /para2/ range,height,rsphfld,rsphcham
      common /para3/ nzeta,npsi,zetmin,zetmax,psimin,psimax
                                     
      ang=angle*pi/180.                 
      dz = (zetmax-zetmin)/(nzeta-1)
      ds = (psimax-psimin)/(npsi-1)
      
C
C  ITERATE TO CALCULATE EACH ELEMENT OF B MATRIX.
C          
      do 5 id =1 1,ndata
      do 5 i=1,16
      do 5 j=1,16
   5     sumdsq(i,j,id) = (0.0,0.0)
            
      do 10 id = 1,ndata                          

      do 10 i=1,np
      do 10 j=1,np    
         ip = 4*(i-1) +j
      do 10 l=1,np
      do 10 m=1,np    
         iq = 4*(l-1) +m           
C
C  DOUBLE INTEGRATION FOR AZ, EL, ANGLES
C               
        sumz = (0.0, 0.0)
        do 20 iz=1,nzeta
          zeta = (iz-1)*dz +zetmin 

          sums = (0.0, 0.0)
        do 30 is=1,npsi
          psi = (is-1)*ds +psimin
                   
C
C  WITH CORRELATION FUNCTIONS 
C
c        sumz1 = (0.0, 0.0)                  { for correlated case below...
c        do 40 jz=1,nzeta
c          zetap = ( (jz-1)*dz +zetmin)
c          sums1 = (0.0, 0.0)
c        do 50 js=1,npsi
c          psip = ( (js-1)*ds +psimin)
c          xxp = tan(psi)/cos(ang+zeta)-tan(psip)/cos(ang+zetap) 
c          yyp = tan(ang+zeta) - tan(ang+zetap)
c          rho = height * sqrt(xxp**2 + yyp**2)
c          gamma = exp(-rho/correl)
c          arg1 = 1. /cos(ang+zeta) /cos(psi)
c          arg2 = 1. /cos(ang+zetap) /cos(psip)
c          z1 = cexp(cj *2.0 *height *q0 *(arg1-arg2))
c          denom = cos(ang+zeta) * cos(ang+zetap)         
c          dsq = d(iz,is,id,i,l) * conjg( d(jz,js,id,j,m))
c          funb = dsq *gamma *z1 /denom   
c        sums1 = sums1 + funb * ds 
c   50   continue                  
c        sumz1 = sumz1 + sums1 * dz
c   40   continue 
C
C  WITHOUT CORRELATION FUNCTION
C                    
c     for az. over el.
c          trig = cos(ang+zeta) * cos(psi)**2
c     for el. over az.
          trig = cos(zeta)**2 
     &               *abs(cos(ang)*cos(psi) -sin(ang)*tan(zeta))

          dsq = d(iz,is,id,i,l) * conjg( d(iz,is,id,j,m))
                                       
          funb = dsq *trig 

          sumdsq(ip,iq,id) = sumdsq(ip,iq,id) +dsq
          sums = sums + funb * ds
   30   continue

          sumz = sumz + sums * dz
   20   continue

        b(ip,iq,id) = sumz *rsphfld**4 /height**2
   10 continue                                    
c      do 80 id=1,ndata
c         print*,(cabs(sumdsq(i,i,id)),i=1,16)
c   80 continue
c       print *,'rsphfld=',rsphfld

      return
      end                                                                 
C                                                                         
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine yvectorc(it)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                                                  
      complex b(16,16,21),sums,sumz,funb
      complex esoil(100,4,21),y(16,21),yy(16,21)
      complex z1,z2,cj 
      real sig(4,21),mueller(4,4),lambda(4,4)
                           
      common /target/ esoil
      common /measur/ y,yy
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip

C
C  CALCULATE AVERAGE VALUES OF MEASURED SURFACE BACKSCATTER INTENSITIES.
C                     
      nm = np**2
      do 10 i=1,nm
      do 10 id=1,ndata                                      
         y(i,id)=(0.0,0.0)
   10 continue
C
      do 20 i = 1,np
      do 20 j = 1,np
         ip = 4*(i-1) +j
      do 20 id = 1,ndata

          y(ip,id) = y(ip,id) + esoil(it,i,id)
     &              * conjg(esoil(it,j,id))   

   20 continue      
        do 30 id=1,ndata
        do 30 i=1,nm
          yy(i,id) = yy(i,id) + y(i,id)/ntrace
   30   continue    
      return
      end
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine solveqc
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                                                  
      complex b(16,16,21),y(16,21),x(16,21),yy(16,21)
      complex z1,z2,cj        
      complex b1(16,16),y1(16,1),x1(16,1)
      
      common /calmat/ b  
      common /measur/ y,yy
      common /xcorr/ x
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
         
      nm=16
C
C  ITERATE FOR EACH FREQUENCY
C                                    
      do 10 id=1,ndata     

        do 20 i=1,nm
          y1(i,1) = y(i,id)
        do 20 j=1,nm
          b1(i,j) = b(i,j,id)
   20   continue      
       
C
C  INVERSE MATRIX B1
C                   
        call matinvc(b1,nm)                             
C
C  MULTIPLY MATRICES B1*Y1=X1
C
        call matmltc(b1,y1,x1,nm,nm,1)                     
                         
C
C AVERAGING W MATRIX
C
        do 30 i=1,nm
c       x(i,id) = x(i,id) + x1(i,1)/ntrace  
	x(i,id) = x1(i,1)
   30   continue             
   10 continue
      return
      end
C
C+++++++++++++++++++++++++++++++++++++++++++
C      subroutine matinvc
C+++++++++++++++++++++++++++++++++++++++++++
C
CI  INVERSION OF COMPLEX MATRIX C, OF ORDER OF LL.
CI INVERSE IS RETURNED IN PLACE OF C.
      SUBROUTINE MATINVC(C,LL)
      COMPLEX C(1),STOR,STO,ST,S
      DIMENSION LR(77)
      COMPLEX X
      CABQ(X)=REAL(X)*REAL(X)+AIMAG(X)*AIMAG(X)
      DO 20 I=1,LL
       LR(I)=I
   20 CONTINUE
      M1=0
      DO 18 M=1,LL
       K=M
       DO 2 I=M,LL
        K1=M1+I
        K2=M1+K
        IF(CABQ(C(K1))-CABQ(C(K2))) 2,2,6
   6   K=I
   2   CONTINUE      
       LS=LR(M)
       LR(M)=LR(K)
       LR(K)=LS
       K2=M1+K
       STOR=C(K2)
       J1=0
       DO 7 J=1,LL
        K1=J1+K
        K2=J1+M
        STO=C(K1)
        C(K1)=C(K2)
        C(K2)=STO/STOR
        J1=J1+LL
   7   CONTINUE
       K1=M1+M
       C(K1)=1./STOR
       DO 11 I=1,LL
        IF(I-M) 12,11,12
   12   K1=M1+I
        ST=C(K1)
        C(K1)=(0.,0.)
        J1=0
        DO 10 J=1,LL
         K1=J1+I
         K2=J1+M
         C(K1)=C(K1)-C(K2)*ST
         J1=J1+LL
   10   CONTINUE
   11  CONTINUE
       M1=M1+LL
   18 CONTINUE
      J1=0
      DO 9 J=1,LL
       IF(J-LR(J)) 14,8,14
   14  LRJ=LR(J)
       J2=(LRJ-1)*LL
   21  DO 13 I=1,LL
        K2=J2+I
        K1=J1+I
        S=C(K2)
        C(K2)=C(K1)
        C(K1)=S
   13  CONTINUE
       LR(J)=LR(LRJ)
       LR(LRJ)=LRJ
       IF(J-LR(J)) 14,8,14
   8   J1=J1+LL
   9  CONTINUE
      RETURN
      END
C
C+++++++++++++++++++++++++++++++++++++++++++
C      subroutine matmltc
C+++++++++++++++++++++++++++++++++++++++++++
CI C=A*B:L*N,  A:L*M,  B:M*N
C
      SUBROUTINE MATMLTC(A,B,C,L,M,N)
      COMPLEX A(L,M),B(M,N),C(L,N)
      DO 10 I=1,L
      DO 10 J=1,N
   10   C(I,J)=(0.,0.)
      DO 15,I=1,L
      DO 15 J=1,N  
      DO 15 K=1,M
        C(I,J)=C(I,J)+A(I,K)*B(K,J)
   15 CONTINUE
      RETURN
      END                                       
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine sigmac
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      complex x(16,21),mueller(4,4),mu(4,4)
      complex z1,z2,cj,muel(21,4,4)
      complex nu(4,4),nuinv(4,4),w(4,4),w1(4,4)
      real sig(4,21),sigav(4),alpha,eta,lambda(4,4)
      real alphax,etax,lambdax(4,4)
      real alp(21),zet(21),alpx(21),zetx(21)
	real max1,max2,min1,min2
	integer outliers(4,4)
      
      common /xcorr/ x   
      common /muell/ mueller,muel
      common /scatt/ sig,sigav
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /phaze/ alpha,eta,alphax,etax
      common /fas/ alp,zet,alpx,zetx
	common /badpts/ outliers

C
C  BACKSCATTERING COEFFICIENTS FOR VV,VH,HV,HH POL.
C
      do 10 id=1,ndata
        sig(1,id) = cabs(x(1,id))
        sig(2,id) = cabs(x(11,id))
        sig(3,id) = cabs(x(6,id))
        sig(4,id) = cabs(x(16,id))
   10 continue    
c
c  AVERAGED BACKSCATTERING COEFFICIENTS
c   (Skip nskip points on each end plus do not include 2 highest
c    and two lowest points)
c           
      do 15 ip=1,4 
	max1 = 0.0
	max2 = 0.0
	min1 = 10.0
	min2 = 10.0
        sigav(ip)=0.0
	do 11, id = 1+nskip,ndata-nskip
		if(sig(ip,id).lt.min1)then
			min1=sig(ip,id)
			outliers(ip,1) = id
		endif
		if(sig(ip,id).gt.max1)then
			max1=sig(ip,id)
			outliers(ip,3)=id
		endif
11	continue
	do 12, id = 1+nskip,ndata-nskip	
		if(id.ne.outliers(ip,1).and.sig(ip,id).lt.min2)then
			min2=sig(ip,id)
			outliers(ip,2)=id
		endif
		if(id.ne.outliers(ip,3).and.sig(ip,id).gt.max2)then
			max2=sig(ip,id)
			outliers(ip,4)=id
		endif
12	continue
	do 13, id = 1+nskip,ndata-nskip	
		if(id.ne.outliers(ip,1).and.id.ne.outliers(ip,2))then
			if(id.ne.outliers(ip,3).and.id.ne.outliers(ip,4))then
			sigav(ip) = sigav(ip) + sig(ip,id)
			endif
		endif
13	continue
	sigav(ip) = sigav(ip)/(ndata - 2*nskip -4)
15	continue

c
c  CONVERT TO DECIBEL UNIT
c
      do 17 ip=1,4                              
        sigav(ip) = 10.*alog10(4.*pi*sigav(ip))
      do 17 id=1,ndata
        sig(ip,id) = 10.*alog10(4.*pi*sig(ip,id))
   17 continue    
C
C  DEFINE NU(INV), NU(INV) VECTORS
C
      do 20 i=1,4
      do 20 j=1,4
        nu(i,j) = (0.0, 0.0)   
        nuinv(i,j) = (0.0, 0.0)   
        mueller(i,j) = (0.0, 0.0)
      do 20 id=1,ndata
        muel(id,i,j)=(0.0,0.0)
   20 continue               

      nu(1,1) = (1.0, 0.0)     
      nu(2,2) = (1.0, 0.0)     
      nu(3,3) = (1.0, 0.0)     
      nu(3,4) = (1.0, 0.0)     
      nu(4,3) = (0.0,-1.0)     
      nu(4,4) = (0.0, 1.0)           
                                   
      nuinv(1,1) = (1.0, 0.0)     
      nuinv(2,2) = (1.0, 0.0)     
      nuinv(3,3) = (0.5, 0.0)     
      nuinv(4,3) = (0.5, 0.0)     
      nuinv(3,4) = (0.0, 0.5)     
      nuinv(4,4) = (0.0,-0.5)     

C
C  GENERATE MODIFIED STOKES SCATTERING OPERATER FROM X VECTOR
C                      
       ntotal=ndata -2*nskip
      do 30 id=1+nskip, ndata-nskip
                       
        w(1,1) = x(1,id)
        w(1,2) = x(6,id)
        w(1,3) = x(2,id)
        w(1,4) = x(5,id)
        w(2,1) = x(11,id)
        w(2,2) = x(16,id)
        w(2,3) = x(12,id)
        w(2,4) = x(15,id)
        w(3,1) = x(3,id)
        w(3,2) = x(8,id)
        w(3,3) = x(4,id)
        w(3,4) = x(7,id)
        w(4,1) = x(9,id)
        w(4,2) = x(14,id)
        w(4,3) = x(10,id)
        w(4,4) = x(13,id)
                      
        call matmltc (nu,w,w1,4,4,4)
        call matmltc (w1,nuinv,mu,4,4,4)

        do 40 i=1,4
        do 40 j=1,4
           mueller(i,j) = mueller(i,j) + 4.0 *pi *mu(i,j)/ntotal
           muel(id,i,j) = 4.0 *pi *mu(i,j)
   40   continue
   30 continue
C
C  MUELLER MATRIX CALCULATION FOR PHASE STATISTICS
C
      LAMBDA(1,1)=cabs(MUELLER(1,1)/2.0)
      LAMBDA(3,3)=cabs(MUELLER(2,2)/2.0)
      LAMBDA(1,3)=cabs((MUELLER(3,3)+MUELLER(4,4))/4.0)
      LAMBDA(1,4)=cabs((MUELLER(3,4)-MUELLER(4,3))/4.0)
      ALPHA=(LAMBDA(1,3)**2+LAMBDA(1,4)**2)/(LAMBDA(1,1)*LAMBDA(3,3))
      ALPHA=SQRT(ALPHA)
      ETA=ATAN(LAMBDA(1,4)/LAMBDA(1,3)) *180./PI
 
c      print*, 'Alpha=',ALPHA,'   Eta=',ETA
C  
      LAMBDAX(1,1)=cabs(MUELLER(1,1))/2.0
      LAMBDAX(3,3)=cabs(MUELLER(1,2))/2.0
      LAMBDAX(1,3)=cabs(MUELLER(1,3))/2.0
      LAMBDAX(1,4)=cabs(MUELLER(1,4))/2.0
      ALPHAX=(LAMBDAX(1,3)**2+LAMBDAX(1,4)**2)
     &                       /(LAMBDAX(1,1)*LAMBDAX(3,3))
      ALPHAX=SQRT(ALPHAX)
      ETAX=ATAN(LAMBDAX(1,4)/LAMBDAX(1,3)) *180./PI

c      print*, 'Alphax=',ALPHAx,'   Etax=',ETAx
C
C  PHASE STATISTICS FOR EACH FREQUENCY
C                                     
      DO ID=1+nskip,NDATA-nskip
        LAMBDA(1,1)=cabs(MUEL(ID,1,1)/2.0)
        LAMBDA(3,3)=cabs(MUEL(ID,2,2)/2.0)
        LAMBDA(1,3)=cabs((MUEL(ID,3,3)+MUEL(ID,4,4))/4.0)
        LAMBDA(1,4)=cabs((MUEL(ID,3,4)-MUEL(ID,4,3))/4.0)
        ALP1=(LAMBDA(1,3)**2+LAMBDA(1,4)**2)/(LAMBDA(1,1)*LAMBDA(3,3))
        ALP(ID)=SQRT(ALP1)
        ZET(ID)=ATAN(LAMBDA(1,4)/LAMBDA(1,3)) *180./PI
 
C                                             
        LAMBDAX(1,1)=cabs(MUEL(ID,1,1))/2.0
        LAMBDAX(3,3)=cabs(MUEL(ID,1,2))/2.0
        LAMBDAX(1,3)=cabs(MUEL(ID,1,3))/2.0
        LAMBDAX(1,4)=cabs(MUEL(ID,1,4))/2.0
        ALPX1=(LAMBDAX(1,3)**2+LAMBDAX(1,4)**2)
     &                       /(LAMBDAX(1,1)*LAMBDAX(3,3))
        ALPX(ID)=SQRT(ALPX1)
        ZETX(ID)=ATAN(LAMBDAX(1,4)/LAMBDAX(1,3)) *180./PI
c      print *,'id,a,z=',id,alp(id),zet(id),alpx(id),zetx(id)

      ENDDO
      return
      end
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      subroutine prtoutc
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                                                  
      complex s0(21),esph0(4,21),esphd(15,15,4,21)
      complex esoil(100,4,21),cj
      complex b(16,16,21),y(16,21),yy(16,21),yav(16)
      complex x(16,21),xav(16)
      complex d(15,15,21,4,4),mueller(4,4),muel(21,4,4)
      complex frt,falpha,fbeta,facrt(21),facalp(21),facbet(21)
      real sig(4,21),sigav(4) 
      real alp(21),zet(21),alpx(21),zetx(21)
      character*30 msoutfile,fasoutfile,gndfile,reffile,smeafile
	integer outliers(4,4)
	character*7,pwords(4)

      common /sphere/ s0,esph0,esphd
      common /target/ esoil
      common /distor/ d  
      common /calmat/ b
      common /measur/ y,yy
      common /xcorr/ x  
      common /muell/ mueller,muel   
      common /scatt/ sig,sigav
      common /phaze/ alpha,eta,alphax,etax   
      common /fas/ alp,zet,alpx,zetx
      common /para1/ ntrace,ndata,np,angle,fstart,fstop,cj,pi,nskip
      common /para2/ range,height,rsphfld,rsphcham
      common /para3/ nzeta,npsi,zetmin,zetmax,psimin,psimax
      common /file1/ msoutfile,fasoutfile,gndfile,reffile,smeafile,it
      common /badpts/ outliers 
                  
C
C  PRINT OUT 
C              
c        write(20,*) 'Sphere (Field) filename= ',reffile
c        write(20,*) 'Sphere (Dist.) filename= ',smeafile
c        write(20,*) 'Soil Backscatter filename= ',gndfile

c        write(21,*) 'Sphere (Field) filename= ',reffile
c        write(21,*) 'Sphere (Dist.) filename= ',smeafile
c        write(21,*) 'Soil Backscatter filename= ',gndfile
c                   
c        write(20,300)     
c        df=(fstop-fstart)/(ndata-1)
c      do 20 id=1,ndata
c        write(20,100) fstart+(id-1)*df,sig(1,id),sig(4,id)
c     &                             ,sig(2,id),sig(3,id)
c   20 continue     
c        write(20,310) nskip 

	if(it .eq. 1) then
	write(20,*) ntrace 
	else
	continue
	endif

        write(20,110) sigav(1),sigav(4),sigav(2),sigav(3)
c	pwords(1) ='For VV:'
c	pwords(2) = 'For VH:'
c	pwords(3) = 'For HV:'
c	pwords(4) = 'For HH:'
c	write(20,311)
c	do 17, i = 1,4
c		write(20,312)pwords(i),(fstart+(outliers(i,j)-1)*df,
c     #			j = 1,4) 
c17	continue
                                          
c        write(20,410)
c      do 56 i =1, 16 
c        write(20,415) i, x(i,6),cabs(x(i,6)),phasec(x(i,6))
c   56 continue       
c          
c      do 25 i=1,16      
c        yav(i)=0.0
c        xav(i) = 0.0
c      do 25 id=1+nskip,ndata-nskip
c         yav(i)=yav(i)+yy(i,id)/(ndata-2*nskip)
c         xav(i)=xav(i)+x(i,id)/(ndata-2*nskip)
c   25 continue  
                     
c        write(20,411)
c      do 66 i =1, 16 
c        write(20,415) i, xav(i),cabs(xav(i)),phasec(xav(i))
c   66 continue       

c        write(21,140)     
        write(21,150) alpha,eta,alphax,etax
c
c  AVERAGED ALPHA AND ZETA
c            
      alpav=0.0
      zetav=0.0
      alpxav=0.0
      zetxav=0.0
      do 15 id=1+nskip,ndata-nskip
         alpav=alpav+alp(id)/(ndata-2*nskip)
         zetav=zetav+zet(id)/(ndata-2*nskip)
         alpxav=alpxav+alpx(id)/(ndata-2*nskip)
         zetxav=zetxav+zetx(id)/(ndata-2*nskip)
   15 continue       
c
c        write(21,120)
c      do 10 id=1+nskip,ndata-nskip
c         write(21,160) fstart+(id-1)*df,alp(id),zet(id)
c     &             ,alpx(id),zetx(id)
c  10 continue
c      write(21,130) alpav,zetav,alpxav,zetxav

c      do 16 id=1+nskip,ndata-nskip,7
c        write(21,205)  fstart+(id-1)*df
c      do 16 i=1,np
c        write(21,260)  i,(cabs(muel(id,i,j)),j=1,np)
c   16 continue          

c  100 format(f7.3,4(f9.3))
  110 format(4(f8.3))   
c  120 format('Freq #',2x,'Alpha,eta for Co-pol',3x,
c     &          'Alpha,eta for X-pol')              
c  130 format('Averaged values',/,2(2x,f8.4,f8.2))
c  130 format(2(2x,f8.4,f8.2))
c  140 format('Alpha,eta for Co-pol',5x,
c     & 'Alpha,eta for X-pol. (from Averaged Mueller)')
  150 format(2(2x,f8.4,f8.2))
c  160 format(f6.3,2(2x,f8.4,f8.2))
c  205 format(/10('*'),' Mueller - matrix(',f5.2,' GHz) ',10('*'))
c  260 format(i4,4f10.6)
c  300 format(/10('*'),' Sigma (dB) ',10('*')/
c     &       '  freq,   VV   -   HH   -   VH   -   HV')    
c  310 format(/'Averaged sigmac (dB) w.r.t. frequency (skip:',i3,')'/)
c  311 format(/'Omitted in averaging following F-pts:')
c  312 format(3x,a7,4(2x,f7.3))
c  410 format(/10('*'),'  X  ',10('*')/                     
c     &       '  No,  Real, Imag., mag., phasec at 5.3 GHz')    
c  411 format(/10('*'),'  X(averaged)  ',10('*')/                     
c     &       '  No,  Real, Imag., mag., phasec at 5.3 GHz')    
c  415 format(i4,3x,2e12.4,3x,f12.5,f10.2)

      return
      end
C+++++++++++++++++++++++++++++++++The end.