        subroutine subcalx(ntrace,type,ang,label,title)
c
C   THIS PROGRAM CALIBRATES THE SCATTEROMETERS USING THE NEW METHOD. 
C
        COMPLEX C2,S11AV,S22AV,S12AV,S21AV,A
        DIMENSION S11AV(21),S22AV(21),S12AV(21),S21AV(21)
        COMPLEX EVV0,EHH0,EVH0,EHV0
        DIMENSION EVV0(21),EHH0(21),EVH0(21),EHV0(21)
C
        COMPLEX EVVS,EHHS,EVHS,EHVS
        DIMENSION EVVS(21),EHHS(21),EVHS(21),EHVS(21)
        COMPLEX EVVU,EHHU,EVHU,EHVU  
        DIMENSION EVVU (120,21),EHHU (120,21)
        DIMENSION EVHU (120,21),EHVU (120,21)  
        COMPLEX S11,S22,S12,S21
        DIMENSION S11 (120,21),S22 (120,21)
        DIMENSION S12 (120,21),S21 (120,21)
C
        CHARACTER*80 REFFILE,SCALFILE,SKYFILE,GNDFILE
        CHARACTER*80 msoutfile 
	character type*1,ang*2,label*11,title*4

        REAL SIG11,SIG22,SIG12,SIG21
        DIMENSION SIG11(21),SIG22(21),SIG12(21),SIG21(21)
        REAL STS11,STS22,STS12,STS21
        DIMENSION STS11(21),STS22(21),STS12(21),STS21(21)
        COMPLEX S0(21),FACT,NORM,z1,z2
        INTEGER INDEX(120),j1,j2   
        
           
      scalfile='~/wheat.dir/data.dir/util.lib/sph_14_x'
      reffile=title//'_ref.'//type
      skyfile =title//'_'//ang//'s.'//type
      gndfile =type//label       
      msoutfile='ax'//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=9.30
      fstop=9.70
      ntrace=ntrace
      npts=21
      ndata=21
      nskip=2

        open(9,file=SCALFILE)
        open(10,file=REFFILE)
        open(12,file=SKYFILE) 
        open(13,file=GNDFILE) 
        PI=4.*ATAN(1.)   
C        ndata=21     
C        ntrace=99
C   READ THEORETICAL VALUE OF SPHERE
        DO 700 I=1,21
        READ(9,*)DUM,S0(I)
700     CONTINUE
C  SPHERE DATA         
        READ(10,103)
         DO 701 I=1,ndata
         READ(10,*)evh0(I)
701     CONTINUE
C
        READ(10,103)
        DO 702 I=1,ndata
        READ(10,*)evv0(I)
702     CONTINUE
        READ(10,103)
         DO 703 I=1,ndata        
         READ(10,*)ehh0(I)
703     CONTINUE
        READ(10,103)
        DO 704 I=1,ndata
        READ(10,*)ehv0(I)
704     CONTINUE
c        do i=1,ndata
c        print *,20*alog10(cabs(evv0(i))),20*alog10(cabs(ehv0(i))),
c     &          20*alog10(cabs(evh0(i))),20*alog10(cabs(ehh0(i))) 
c        enddo
C SKY DATA   
C
        READ(12,103)
        DO 705 I=1,ndata
         READ(12,*)evhs(i)
705     CONTINUE
        READ(12,103)
        DO 706 I=1,ndata
        READ(12,*)evvs(i)
706     CONTINUE
        READ(12,103)
         DO 707 I=1,ndata
         READ(12,*)ehhs(i)
707     CONTINUE
        READ(12,103)
        DO 708 I=1,ndata
        READ(12,*)ehvs(i)
708     CONTINUE
c
c  Set sky to zero; ie, don't use sky subtraction.
c
c        do 832, I = 1,ndata
c        evhs(i) = (0.,0.)
c        evvs(i) = (0.,0.)
c        ehhs(i) = (0.,0.)
c        ehvs(i) = (0.,0.)
c832     continue
c
C UNKNOWN DATA
        DO 722 I=1,ndata
        SIG11(I)=0.
        SIG22(I)=0.
        SIG12(I)=0.
        SIG21(I)=0.
        STS11(I)=0.
        STS22(I)=0.
        STS12(I)=0.
        STS21(I)=0.
722     CONTINUE
      DO 713 ITRACE=1,ntrace
        READ(13,103)
103     FORMAT(1X)
         DO 709 I=1,ndata
         READ(13,*)EVHU(ITRACE,I)
         EVHU(ITRACE,I)=EVHU(ITRACE,I)-EVHS(I)
709     CONTINUE
        READ(13,103)
        DO 710 I=1,ndata
        READ(13,*)EVVU(ITRACE,I)
        EVVU(ITRACE,I)=EVVU(ITRACE,I)-EVVS(I)
710     CONTINUE
        READ(13,103)
         DO 711 I=1,ndata
         READ(13,*)EHHU(ITRACE,I)
         EHHU(ITRACE,I)=EHHU(ITRACE,I)-EHHS(I)
711     CONTINUE
        READ(13,103)
        DO 712 I=1,ndata
        READ(13,*)EHVU(ITRACE,I)
        EHVU(ITRACE,I)=EHVU(ITRACE,I)-EHVS(I)
712     CONTINUE
713     CONTINUE
c
       close(9)
       close(10)
       close(12)
       close(13)
C        DF=(1.4-1.1)/(ndata-1)   
        DF = (FSTOP-FSTART)/(NDATA-1)
C this is a calibration based on c1=c2
        DO 714 ITRACE=1,NTRACE
        DO 715 I=1,ndata
        A=EVH0(I)*EHV0(I)/(EVV0(I)*EHH0(I))
        S12(ITRACE,I)=0.5*S0(I)*(EVVU(ITRACE,I)
     &          /EVV0(I)+EHHU(ITRACE,I)/EHH0(I)
     &          -2.*EVHU(ITRACE,I)/EVH0(I))*(CSQRT(A)/(A-1))
        S21(ITRACE,I)=0.5*S0(I)*(EVVU(ITRACE,I)
     &         /EVV0(I)+EHHU(ITRACE,I)/EHH0(I)
     &         -2.*EHVU(ITRACE,I)/EHV0(I))*(CSQRT(A)/(A-1))
        S11(ITRACE,I)=0.5*S0(I)*( (1/(1-A)+1./CSQRT(1.-A))*
     &      EVVU(ITRACE,I)/EVV0(I)
     &      +(1/(1-A)-1./CSQRT(1.-A))*EHHU(ITRACE,I)/EHH0(I)
     &      +2.*A*EVHU(ITRACE,I)/(EVH0(I)*(A-1)) )
        S22(ITRACE,I)=0.5*S0(I)*( (1/(1-A)+1./CSQRT(1.-A))*
     &       EHHU(ITRACE,I)/EHH0(I)
     &      +(1/(1-A)-1./CSQRT(1.-A))*EVVU(ITRACE,I)/EVV0(I)
     &      +2.*A*EVHU(ITRACE,I)/(EVH0(I)*(A-1)) )  
715     CONTINUE
C        PRINT *,ITRACE 
714     CONTINUE
        RANGE = 12.0 - 1.15
C
C  Following section selects a value or 'factor' to be used in the
c  illumination integral.  These values were computed considering the
c  antenna pattern out to the -15dB level.
c
        IF(ANGLE.EQ.20.0)THEN
                FACT_0 = 1.0/(.66E-4*RANGE**4)
        ELSEIF(ANGLE.EQ.30.0)THEN
                FACT_0 = 1.0/(.71E-4*RANGE**4)
        ELSEIF(ANGLE.EQ.40.0)THEN
                FACT_0 = 1.0/(.808E-4*RANGE**4)
        ELSEIF(ANGLE.EQ.50.0)THEN
                FACT_0 = 1.0/(.955E-4*RANGE**4)
        ELSEIF(ANGLE.EQ.60.0)THEN
                FACT_0 = 1.0/(.125E-3*RANGE**4)     
        ELSEIF(ANGLE.EQ.70.0)THEN
                FACT_0 = 1.0/(.185E-3*RANGE**4)

        ENDIF
c
        open(15,file=msoutfile)
        DO 716 I=1,NDATA    
        DO 717 ITRACE=1,NTRACE
        SIG11(I)=SIG11(I)+CABS(S11(ITRACE,I))**2/ntrace
        SIG22(I)=SIG22(I)+CABS(S22(ITRACE,I))**2/ntrace
        SIG12(I)=SIG12(I)+CABS(S12(ITRACE,I))**2/ntrace
        SIG21(I)=SIG21(I)+CABS(S21(ITRACE,I))**2/ntrace
717     CONTINUE
716     CONTINUE
        slav11=0.
        slav22=0.
        slav12=0.
        slav21=0.                           
C
C  SKIP THE FREQUENCY POINTS OF BOTH ENDS
C
C        NSKIP=3
        NDT=NDATA-(2*NSKIP)
        DO 718 I=1+NSKIP,NDATA-NSKIP
         slav11=slav11+sig11(i)/NDT
         slav22=slav22+sig22(i)/NDT
         slav12=slav12+sig12(i)/NDT
         slav21=slav21+sig21(i)/NDT
718     CONTINUE 
	 write(15,*) ntrace
        DO 719 I=1,ndata
         F=(I-1)*DF+1.1
         WRITE(15,810)F,10.*ALOG10(4.*PI*FACT_0*SIG11(I))
     &               ,10.*ALOG10(4.*PI*FACT_0*SIG22(I))   
     &               ,10.*ALOG10(4.*PI*FACT_0*SIG12(I))
     &               ,10.*ALOG10(4.*PI*FACT_0*SIG21(I))   
719     CONTINUE
810     format(1X,f5.3,5x,4(f8.4,3x))
         WRITE(15,811)'Average=',10.*ALOG10(4.*PI*FACT_0*SLAV11)
     &               ,10.*ALOG10(4.*PI*FACT_0*SLAV22)   
     &               ,10.*ALOG10(4.*PI*FACT_0*SLAV12)
     &               ,10.*ALOG10(4.*PI*FACT_0*SLAV21)   
811     format(1X,A8,5x,4(f8.4,3x))
        return
        END
C*******************************************************************
C
        FUNCTION PHASEX(Z)
C
C******************************************************************
        COMPLEX Z 
        PI=4.*ATAN(1.)
        X=REAL(Z)
        Y=AIMAG(Z)
        PHASE=(180./PI)*ATAN2(Y,X)
        RETURN
        END

