program extract

c  compile with ...
c  g77 -ffixed-line-length-none -g -o Extract_RD_Stats  Extract_RD_Stats.f 

      implicit none

c  for extracting the pair correlation along the line Y=0

      integer i,j,k,l,m,size
      parameter(size = 1001)
      real array(size,size)
      character filename*50
      integer xsize,ysize
      real rtmp1,rtmp2,sum
      real angle(3)
      real scint(0:6),scleg(0:6),scleg2(0:6),cval(0:10)
      character outfile*60
      integer nval,nbelow0,nscale
      real colrval(3,0:10)
      real x_plotsize,y_plotsize
      real ave
      real wdth,rmax
      REAL SMF((-1*size):size)
      real red,green,blue
      integer exit_status
      integer dot
      logical done

      integer ispec,ioffpp,spvall,ilegg,ilabb,nhii,ndeccn,nlbll
      integer mscall,ldsh,hgtlab
      common/conpar/ispec,ioffpp,spvall,ilegg,ilabb,nhii,ndeccn,nlbll,
     +mscall,ldsh,hgtlab

      integer n,lscal,ldash,ndecn,dimen
      parameter (dimen=200)
      real xlen,ylen
      real testarray(dimen,dimen)

c  CODE:

c****  test section
c$$$      ispec=1
c$$$      ilabb = 1
c$$$      call newdev('testarray.ps',1)
c$$$      call psinit(.true.)
c$$$      nhii=-1
c$$$      do j=1,10 
c$$$         do i=1,10 
c$$$            testarray(i,j)=i*j 
c$$$         enddo
c$$$      enddo
c$$$      xlen=5. 
c$$$      ylen=5. 
c$$$      nval=6 
c$$$      call border(xlen, ylen, 1111, 1111, 4, 2, 5, 3) 
c$$$      do n = 1,nval 
c$$$         cval(n)=(n-1)*20 
c$$$      enddo
c$$$      nlbll = 0 
c$$$      mscall = 1                !  different from example 
c$$$      call conrec(testarray, dimen, 10, 10, xlen, ylen, cval, nval) 
c$$$      ldsh = 6                  !  different
c$$$      ndeccn = -1 
c$$$      do n = 1,nval 
c$$$         cval(n)=10.+(n-1)*20 
c$$$      enddo
c$$$      call conrec(testarray, dimen, 10, 10, xlen, ylen, cval, nval) 
c$$$      call plotnd
   
c*****

      ispec=1
      ilegg=0
      ilabb=1
      nhii=-1
      mscall = 1
c turns off legend, labels, hi/lo labels
      scleg(0) = 0.2
      cval(1) = 0.2
      scleg(1) = .5
      cval(2) = .5
      scleg(2) = 0.8
      cval(3) = 0.8
      scleg(3) = 1.2
      cval(4) = 1.2
      scleg(4) = 2.
      cval(5) = 2.
      scleg(5) = 3.
      cval(6) = 3.
      nbelow0=1
      nscale=1
       
      print*,'This extracts the Correlation stats along Y=0'
      print*,'Enter filename to read (e.g. 7075L_2PtStats.txt)'
      read(*,*) filename
      open(unit=11,file=filename, status='old')
      read(11,*) xsize,ysize
      print*,'x,y sizes',xsize,ysize
      xsize = 1 + 2*xsize
      ysize = 1 + 2*ysize
      print*,'actual x,y sizes',xsize,ysize
      x_plotsize = 5.
      y_plotsize = 5.*float(ysize)/float(xsize)
      if(xsize.gt.size) stop 'X too large'
      if(ysize.gt.size) stop 'Y too large'
      read(11,*) rtmp1
      read(11,*) (angle(i),i=1,3)
      print*,'angles: ', (angle(i),i=1,3)

      sum = 0.
      do i = 1,xsize
         do j = 1,ysize
            read(11,*) rtmp1, rtmp2
            if(abs(rtmp2).lt.1.e-5) then
               array(i,j) = 0.
            else
               array(i,j) = rtmp1/rtmp2
            endif
c            print*,'i  j  rtmp1  rtmp2  array()'
c            print*,i,j,rtmp1,rtmp2,array(i,j)
            sum = sum + array(i,j)
         enddo
      enddo
      ave = sum/float(xsize*ysize)
      print*,'Average level = ',ave
      close(11)

      dot = 0
      done = .false.
      do i = 30,1,-1
         if((.not.done).and.filename(i:i).eq.'.') then
            dot = i-1  ! point to last character before dot
            done = .true.
         endif
      end do

      do i = 1,xsize
         do j = 1,ysize
            array(i,j) = array(i,j)/ave
         enddo
      enddo

c$$$      print*,'every 5th value'
c$$$      do i = 1,xsize,5
c$$$         print"(1x,10(1x,f8.3))",(array(i,j),j=1,ysize,5)
c$$$      enddo

c      print*,'Enter filename to write out Y=0 stats'
c      read(*,*) filename
      open(unit=11,file=filename(1:dot)//'_Y0.txt', status='unknown')
      write(11,*) '  X  2Pt_stats'
      do i = 1,xsize
         rtmp1 = (array(i,(ysize/2)) + array(i,(ysize/2)+1)
     1      + array(i,(ysize/2)+2))/3.
         write(11,"(1x,i6,1x,f10.6)") i,rtmp1
      enddo

      close(11)

      if(done) then
         outfile = filename(1:dot)//'_contours.ps'
      else
         print*,'Using default output name 2D-contours.ps'
         outfile='2D-contours.ps'
      endif
      write(*,*) 'Using ',outfile,' as output filename'
      call newdev(outfile,1)
      call psinit(.true.)
      call plot(1.,1.,-3)

c      wdth = 3.0
      write(*,*) 'Enter smoothing width, e.g. 6'
      read(*,*) wdth
      CALL SMARRAY(SMF,WDTH)
      CALL SMoothit(SMF,array,rmax,xsize,ysize)

c     initiating post script plotting
      call border(x_plotsize,y_plotsize,0000,1111,4,9,1,9)


      print*,'every 5th value, after smoothing ...'
      do i = 1,xsize,5
         print"(1x,10(1x,f8.3))",(array(i,j),j=1,ysize,5)
      enddo

      red=0.
      green=0.
      blue = 0.8
      call setcolr(red,green,blue)
      CALL SETLW (.015)
      call conrec(array,size,xsize,ysize
     1    ,x_plotsize,y_plotsize,cval(1),2)
      red = 0.
      green = 0.
      blue = 0.
      call setcolr(red,green,blue)
      CALL SETLW (.015)
      call conrec(array,size,xsize,ysize
     1    ,x_plotsize,y_plotsize,cval(3),2)
      red = 1.
      green=0.
      blue = 0.
      call setcolr(red,green,blue)
      CALL SETLW (.015)
      call conrec(array,size,xsize,ysize
     1    ,x_plotsize,y_plotsize,cval(5),2)
      call plotnd

      if(done) then
         call system('convert '//outfile//' '//filename(1:dot)
     1        //'_contours.jpg',exit_status)
      else
         call system('convert '//outfile//' 2D-contours.jpg',exit_status)
      endif
      if(exit_status.eq.0) then
         if(done) then
            call system('rm '//filename(1:dot)//'_contours.ps')
         else
            call system('rm 2D-contours.ps')
         endif
      endif

      call exit(0)
      end

c
c ________________________________________
c
        SUBROUTINE SMARRAY(SMF,RWDTH)
c  WDTH is the width of the smoothing filter
c  SMF is the array that contains the filter values
c  from JS Kallend's smoothing program for pole figures!
c
c        IMPLICIT INTEGER*2 (I-N)
	implicit none

	integer i,k,nnsize

	parameter (nnsize = 1001)

	REAL SMF((-1*nnsize):nnsize),dif,tot,wdth,rwdth,deg,sdeg

c  CODE::

c        DO 100 I=2,19
c	do 100 i=1,19
c        DEG=(I-1)*0.087266
c        SDEG=SIN(DEG)
C       5*PI/180
c        RWDTH=WDTH/SDEG
c  special version: adr: no adjustment of width for declination angle
c  dec 01
c
c	write(*,*) 'Enter smoothing width, e.g. 5'
c	read(*,*) wdth
c	RWDTH=WDTH
c
c       PRINT *,I,RWDTH,SDEG

        DO 200 K=0,nnsize
c			DIF=K*5./RWDTH
           DIF = float(K)/RWDTH   !  take out 5 degree grid assumption
           IF(DIF.GT.8.)THEN
              SMF(K) = 0.
              SMF(-1*K) = 0.
           ELSE
              SMF(K) = EXP(-1.*DIF*DIF)
              SMF(-1*K) = SMF(K)
           ENDIF
c           print *,'index, interm. result= ',k,smf(k),-1*k,smf(-1*k)
 200    continue
c
        TOT = 0.
        DO 210 K = 1-nnsize,nnsize
           TOT = TOT+SMF(K)
 210    CONTINUE

        print*,'integrated intensity in SMF = ',tot
 
        DO 220 K = -1*nnsize,nnsize
          SMF(K) = SMF(K)/TOT
 220    CONTINUE
 100    CONTINUE
c  now for special treatment of the first row
c        DO 230 K=-1*nnsize,nnsize
c230     SMF(1,K)=1./72
c
c	DO 300 I=1,nnsize
        SMF(-1*nnsize)=SMF(-1*nnsize)/2.
        SMF(nnsize)=SMF(-1*nnsize)

c        print*,'SMF = '
c        PRINT 1000,(SMF(K),K=(-1*nnsize),nnsize)
 1000   FORMAT('smf=',1X, 10F8.3)
c300	CONTINUE

        RETURN
        END

c ________________________________________
c

        SUBROUTINE SMoothit(SMF,array,rmax,xsize,ysize)

        implicit none

        integer nnsize, xsize,ysize
	parameter (nnsize = 1001)

c	IMPLICIT INTEGER*2 (I-N)
c        INTEGER*2 ICOUNT(19,72)

        integer i,j,k,ix,iy,indx
	real array(nnsize,nnsize)

	REAL SMF((-1*nnsize):nnsize),smct(nnsize)
c        REAL SMF(19,-36:36),SMCT(72)
	real rmax,ccnt
	rmax=0.
        
        print*,'entering SMoothit'
        write(*,*) 'ICOUNT:  '

c$$$	write(*,*) 'SMoothit debug: 1st column, 1st index, every 20'
c$$$	do j=1,nnsize,20
c$$$           print*,'last print: row number ',j
c$$$           PRINT 1000,(array(I,j),i=1,nnsize,20)
c$$$	enddo
c$$$

c	DO 100 I=1,nnsize
	DO I = 1,xsize
c           print*,'1st loop, i= ',i,'  of ',nnsize
           if(mod(i,10).eq.0) print*,'1st loop, i= ',i,'  of ',xsize
c           DO 150 K=1,nnsize
           DO K=1,ysize
              CCNT = 0.
              DO 160 IX = (-1*nnsize),nnsize
                 if(smf(ix).gt.0.005)then
                    INDX = K+IX
c                    IF(INDX.LT.1) INDX = INDX+nnsize
                    IF(INDX.LT.1) INDX = INDX+ysize
c                    IF(INDX.GT.nnsize) INDX = INDX-nnsize
                    IF(INDX.GT.ysize) INDX = INDX-ysize
                    CCNT = CCNT+array(I,INDX)*SMF(IX)
                 endif
 160          CONTINUE
              SMCT(K) = CCNT
c 150       CONTINUE
           enddo

c           DO 110 K=1,nnsize
           DO K = 1,ysize
c     array(I,K)=SMCT(K)+.5
              array(I,K) = SMCT(K)
              if(array(i,k).gt.rmax) rmax = array(i,k)
c 110       CONTINUE
           enddo
c 100    CONTINUE
        enddo
c$$$        print*,'rmax = ',rmax
c$$$
c$$$	write(*,*) 'SMoothit debug: 1st column, 2nd index'
c$$$        do j=1,nnsize,10
c$$$           print*,'row number ',j
c$$$           PRINT 1000,(array(i,j),i=1,nnsize,10)
c$$$c 1000      FORMAT(10(1x,f9.3))
c$$$        enddo

c	DO I=1,nnsize
	DO I = 1,ysize
c           print"(' row number  ',i8,'  of ',i8)",i,nnsize
           if(mod(i,10).eq.0) then
              print"(' 2nd loop:  ',i8,'  of ',i8)",i,ysize
           endif
c           DO K = 1,nnsize
           DO K = 1,xsize
              CCNT = 0.
              DO IX = (-1*nnsize),nnsize
                 if(smf(ix).gt.0.005)then
                    INDX = K+IX
c                    IF(INDX.LT.1) INDX = INDX+nnsize
c                    IF(INDX.GT.nnsize) INDX = INDX-nnsize
                    IF(INDX.LT.1) INDX = INDX+xsize
                    IF(INDX.GT.xsize) INDX = INDX-xsize
                    CCNT = CCNT+array(indx,I)*SMF(IX)
                 endif
              enddo
              SMCT(K) = CCNT
           enddo

c           DO K=1,nnsize
           DO K = 1,xsize
c     array(k,i)=SMCT(K)+.5
              array(k,i)=SMCT(K)
              if(array(k,i).gt.rmax) rmax=array(k,i)
           enddo
c     PRINT 1000,(array(I,K),K=1,nnsize)
1000   FORMAT(18(1x,g10.2))
      enddo

c$$$	write(*,*) 'SMoothit debug: 1st column, 1st index'
c$$$	do j=1,nnsize,10
c$$$           print*,'last print: row number ',j
c$$$           PRINT 1000,(array(I,j),i=1,nnsize,10)
c$$$	enddo

      RETURN
      END

      include '/code/paul.lee.codes/psplot.txt'