!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                    Peak.f90  v1.0.1                        !
!                                   Summer 2012                              !
!        note: functionality has been added several times      	             !
!                                                                            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Release Notes:                                                             !                                            !
!           -fixed issues with minimum threshold. Anything below about 100   !
!             is unsuitble                                                   !
!        "         Gauss subtract:                                           !
!                  Says end of file in RTE no idea why                       !   
!        "         Parameters                                                !
!            -Very buggy but not important and bottom of the priority list   !
!        "         measure peak now measures alpha value                     !
!        "         write... writes tons of stuff to a file, tis good         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


PROGRAM MENU
  Integer :: val
  LOGICAL :: C
  INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  REAL, DIMENSION(0:4,0:60) :: TABLER
  C=.TRUE.
  Tablei=0
  tabler=0
  
  DO WHILE (C.EQV..TRUE.)  
     Print *,"|---------------------------------------------------------------------|"
     Print *,"|                     welcome to SAT v 1.0                            |"
     print *,"|            please select an option from the list below              |"
     Print *,"|---------------------------------------------------------------------|"
     Print *,"|      0  -- Find new peaks                                           |"
     print *,"|      1  -- Measure Peaks growth                                     |"
     print *,"|      2  -- Prepare 1D Gaussian Bins                                 |"
     print *,"|      3  -- Prepare peaks for imaging                                |"
     print *,"|      4  -- Subtract Gaussian from Data [run-time-FDNE]              |"
     print *,"|      5  -- Write Table values to File                               |"
     print *,"|      6  -- Import Table values from File                            |"
     print *,"|      7  -- SORT Alpha Table                                         |"
     print *,"|      8  -- Q Space Gaussian-- pretty buggy                          |"
     print *,"|      9  -- PARAMETERS N/A                                           |"
     print *,"|      10 -- exit                                                     |"
     print *,"|_____________________________________________________________________|"
     Read (*,*) val
     
     If (val.eq.0) THEN
        call Peakfinder(TABLER,TABLEI)
           !goto 23
     else If (val.eq.1) THEN
        Call ARPHA(tabler,tablei)
        !   goto 23
     else If (val.eq.2) THEN
        call Gauss(TABLER,TABLEI)
     else If (val.eq.3) THEN
        call time_ev(TABLER,TABLEI)
        !   goto 23!
     else If (val.eq.4) THEN
        call Gauss_subtract(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
        !   goto 23
     else IF(val.eq.5)THEN
        call WRIT(TABLER,TABLEI)
     else If (val.eq.6) THEN
        call IMPORT(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
     else If (val.eq.7) THEN
        call SORT(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
        !   goto 23
     else If (val.eq.8) THEN
       call Q_SPACE(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
        !   goto 23
     else If (val.eq.10) THEN
        C=.FALSE.
        EXIT
     else
        print *, "please select a valid menu option"
        !   goto 23
     endif
     !Kenny Peou is great,
     !Menu would not have happened,
     !without his input...
  END DO
  STOP
END PROGRAM MENU

!Table format
!TABLEI
! Group# | xpeak | ypeak  | xmin |ymin |xmax |ymax| 
!
!TABLER
!phi peak \ rpeak \ alpa \ dphi 
!
!
SUBROUTINE PEAKFINDER(tabler,tablei)

  !IMPLICIT NONE
  REAL, DIMENSION(0:4,0:60) :: tabler
  INTEGER, DIMENSION(0:6,0:60) :: tablei
  REAL, DIMENSION(0:299,0:299) :: s,q,p
  INTEGER :: x, y, x1, y1, x2, y2, xi, yi, xf, yf, xii, yii, xmax0, x0, y0, xmax1, ymax1, xmin1, ymin1, wxpeak, wypeak
  LOGICAL Spot, Spot2, Newspots, Write, Gauss
  INTEGER,SAVE :: xmin,ymin,xmax,ymax, lambda
  REAL :: xpeak, ypeak, I, I2, Imax, rmax, rmin, r, mint, r0

  !just opening the files and assigning them to an array
  OPEN(923,FILE='/home/ramalab/research/dj11/dj11f355.dat',STATUS='OLD',ACTION='READ')
  OPEN(900,FILE='OUTPUT.data')
  READ(923,*)s
  lambda=0
  !defining out coordinate system
  !r0=sqrt(x0*x0+y0*y0)
  x0=142
  y0=169
  rmin=r0+82
  rmax=rmin+48
  mint=100.0
  Spot=.FALSE.
  Spot2=.FALSE.
  ymin=0
  ymax=0
  xmin=0
  xmax=0
  
  !Note this next line may need ot be moved, I am not sure that this is entirely kosher
  q=s
  p=s
!  WHERE(s>rmin.and.s<rmax)
     !easiest way to do this loop hence why we do this in Fortran not C >.>
     !when we find a point in the two rings that meets this condition of a min value 
  do x=0,299
     do y=0,299
        q(x,y)=p(x,y)
        I=q(x,y)
        if (I >= mint)then
        !   PRINT *,"Found one!"
        !   Print *,lambda
           tablei(0,lambda)=lambda
        !   lambda=lambda+1
        !   Print*,I,x,y
           Spot=.TRUE.
           spot2=.TRUE.
           !variables come
           !also variables go
           !the below are saved
           x1=x
           y1=y
           ymin=y
           xmin=x
           ymax=y
           xmax=x

           p=q
           Imax=0.0
          !PRint *," ENTERING ACTUAL FIRST LOOP"
           DO WHILE(Spot.eqv..TRUE.)
              I=p(x0,y0)
              If(I < mint)THEN
                 Spot=.FALSE.
              endif
              y0=y0+1
              y0=ymax
           END DO
           !PRINT*,"ENDING FIRST LOOP"
           ymin1=ymin
           ymax1=ymax
           !print*,"starting loop1"
           Do y0=ymin1,ymax1
              SPOT2=.TRUE.
              x0=x1
              DO WHILE(spot2.eqv..TRUE.)
                 I=p(x0,y0)
                 if(I>Imax)THEN
                    !some Haikus composed for the issue
                    !---Here is the problem
                    !What it is I do not see
                    !tab is off in loop?---
                    !---Line one-thirty five
                    !One thirty six,between them,
                    !Something goes awry.
                    !------------------------new style-----------Actual Rhyme ---- Ideal
                    !when implicit equals some,                      A              A
                    !values undefined,                               B              B
                    !are no problem of any kind.                     B              A
                    !though if we set implicit none,                 C              B
                    !then we find,                                   B              C
                    !programming becomes more fun.                   C              B
                    !implicit defining creates more problems,        D              C
                    !than it illuminates answers,                    E              A
                    !yeah, I am out of ideas, schmoblems...          D              C
                    xpeak=x0
                    ypeak=y0
                    Imax=I
                 endif
                 if(I<mint)THEN
                    Spot2=.FALSE.
                 endif
                 x0=x0+1
                 if(x0>xmax)THEN
                    xmax=x0
                 endif
              end do
              SPOT2=.TRUE.
              x0=x1
              DO WHILE(spot2.eqv..TRUE.)
                 I=p(x0,y0)
                 if(I>Imax) THEN
                    xpeak=x0
                    ypeak=y0
                 end if
                 if (I<mint)THEN
                    Spot2=.FALSE.
                 endif
                 x0=x0-1
                 if(x0<xmin)THEN
                    xmin=x0
                 endif
              end do
           end do
           
           xmin1=xmin
           xmax1=xmax
           
           Do x0=xmin1,xmax1
              Spot2=.TRUE.
              y0=y1
              DO WHILE(spot2.eqv..TRUE.)
                 I2=p(x0,y0)
              if(I2>Imax) THEN
                 xpeak=x0
                 ypeak=y0
              end if
              if(I2<mint)THEN
                 Spot2=.FALSE.
              end if
              y0=y0+1
              if(y0>ymax)THEN
                 ymax=y0
              endif
           end do
           SPOT2=.TRUE.
           y0=y1
           !Print *,"alg in process pease stand by"
           DO WHILE(spot2.eqv..TRUE.)
              I2=p(x0,y0)
              if(I2>Imax) THEN
                 xpeak=x0
                 ypeak=y0
                 Imax=I
              end if
              IF(I2<mint)THEN
                 Spot2=.FALSE.
              endif
              y0=y0-1
              if(y0<ymin)THEN
                 ymin=y0
              endif
              !print *,ymin
           end do
        end do
        
        do x0=xmin,xmax
           do y0=ymin,ymax
              p(x0,y0)=0.0
           end do
        end do
        If(xmin>2.AND.ymin>2)THEN
           xmin=xmin-2
           ymin=ymin-2
        endif
        if(xmax<297.AND.ymax>297)THEN
           xmax=xmax+2
           ymax=ymax+2
        endif
        wxpeak=NINT(xpeak)
        wypeak=NINT(ypeak)
        Write(900,'(6I4)')xmin, ymin, xmax, ymax, wxpeak, wypeak
        !Print*,xpeak,ypeak,xmin,ymin !we get stuck like 1/3 of the way through, why?

        tablei(1,lambda)=xpeak
        tablei(2,lambda)=ypeak
        tablei(3,lambda)=xmin
        tablei(4,lambda)=ymin
        tablei(5,lambda)=xmax
        tablei(6,lambda)=ymax
       ! print *,tablei(0,lambda)
        !Late july, a Haiku on issues unknown
        !---While I debug this,
        !I begin to take notice,
        !How much is broken. ---
        !--- Array 2 9 9
        !Loop goes zero three hundred
        !Doesn't matter at all---
        !---If statement shouldn't end
        !Add else condition to end
        !Doesn't matter at all---
        !---One third the way through
        !The loop fails to continue,
        !Why on Earth one third?---
        !---Every time you ask?
        !That is the sign of a bug
        !does one third matter?---
        !---or coincidence?
        !hard to believe it does not.
        !Hard to connect dots.
        !~~~~~~~~~~~~~~~~~~~~~~~~~~Let's switch styles~~~~~~~~~~~~~~~~~~~~~~~~~~~
        !Looping over variable x,
        !Looping over variable y,
        !as long as I don't care, try exchange,
        !Y before repeats the same this try.
        ! MY STYLE NEEDS MORE WORK THAN THIS CODE
        
        q=p
        lambda=lambda+1
     elseif(I < mint)THEN
        !do a little dance?
     end if
  end do
end do
close(900)
close(923)
Print*,"That should be the last of them skipper!"
RETURN
END SUBROUTINE PEAKFINDER !(TABLER, TABLEI)



SUBROUTINE time_ev(TABLER,TABLEI)
IMPLICIT NONE
REAL,DIMENSION(0:4,0:60) :: TABLER
INTEGER, DIMENSION(0:6,0:60) :: TABLEI

! DO NOT CHANGE THESE CHANGE THE FILES TO FIT THE STRUCTURE DO NOT CHANGE THESE!
 integer :: i,j,l,m,z,nums,spots
 ! Used to put variables in file names
 CHARACTER(len=1) :: cnum ! 1-9
 CHARACTER(len=2) :: cnum2 ! 10-99
 CHARACTER(len=3) :: cnum3 !100-360
 CHARACTER(len=3) :: posx !xvalue initially found at  --\___might want to change
 CHARACTER(len=3) :: posy !y value initially found at --/    to peak center

! Used to mark directories


 CHARACTER(len=33) :: n1 !/home/michael/research/dj11/dj11f+ other
 CHARACTER(len=41) :: n10 !/home/michael/research/dj11/Spot_ev/dj11f+ other
 CHARACTER(len=37) :: n11 !/home/michael/research/dj11/IvT/dj11f+ other
 CHARACTER(len=40) :: n2 !output path same as input+ other
 CHARACTER(len=4) :: n3 !.dat
 CHARACTER(len=1) :: n6 !used to add 0
 CHARACTER(len=1) :: n7 ! used to add 0
 CHARACTER(len=1) :: u1 ! used to add _
 CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added 
 CHARACTER(len=47) :: n5 !outputs for Intensity vs time

!this one is for file format -- didn't use it
 CHARACTER(len=3) :: daftpunk
 

 REAL, DIMENSION(0:299,0:299) :: S,g
 INTEGER :: x,y,XMAX,YMAX,XMIN,YMIN,XII,YII,x0,y0,dim1,dim2,k,xpeak,ypeak
 REAL ::  Sum,I2,mint
 

!OPENS every file, gets info for 
OPEN(unit=900,FILE='OUTPUT.data',STATUS='old',ACTION='read')

Print *,"Please indicate the number of spots in the file."
Read (*,*) nums

Do spots=1,nums
   Read(900,'(6I3)')xmin,ymin,xmax,ymax,xpeak,ypeak
   ! Stuff I probably don't need
   
   
   X=XMIN
   Y=YMIN
   sum=0
   mint=100.0
   ! Stuff I may need? 
   
   dim1=xmax-xmin
   dim2=ymax-ymin
   
   ! REAL, DIMENSION(0:dim1,0:dim2) :: V
   WRITE(posx,'(i3)')xmin
   WRITE(posy,'(i3)')ymin
   WRITE(daftpunk,'(i3)')dim1 !totally not used
   
   !Definted for all Three Loops
   N1='/home/ramalab/research/dj11/dj11f'
   N10='/home/rabalab/research/dj11/Spot_ev/dj11f'
   N11='/home/ramalab/research/dj11/Ivt/dj11f'
   n3='.dat'
   u1='_'
   
   n5=n11//posx//posy//n3  !37+10=47
   
   m=0
   l=0
   
   
   WRITE(n6,'(i1)')m
   WRITE(n7,'(i1)')l
   
   !I'm just baffled how this won't worked....
   OPEN(23,file=n5,STATUS='NEW',ACTION='WRITE')
   PRINT *,'You are now in the subroutine, enjoy it while it lasts >.>'
   print*,"dim1, dim2, xmin, ymin, xmax, and ymax: "
   print*,dim1,dim2,xmin,ymin,xmax,ymax
   
   
   !------------------------------------~~----------------------------------------!
   !                  These are the loops That Try mens souls...                  !
   !                                                                              !
   !                        First Iteration of 'THE' Loop                         !
   !------------------------------------------------------------------------------!
   
   LOOP1: DO i=1,9
      j=i+400
      
      !ALL THOSE STRINGS AND STUFF ARE A SIMPLE WAY TO OPEN THE FILES THAT WE WANT OPEN. I WILL IMPLIMENT THIS EARLIER IN THE CODE AS WELL
      
      sum=0
      WRITE(CNUM,'(i1)')i
      n2=n1//n6//n7//cnum//n3 !33+7 =40 
      n4=n10//posx//posy//u1//n6//n7//CNUM//n3  !41+6+4+3=55
      !tricky way to make sure we dont have two files open with the same tag
      
      
      open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
      Open(unit=j,file=n4,STATUS='NEW',ACTION='WRITE')
      
      READ(i,*)s
      k=0
      do y=ymin, ymax
         do x=xmin, xmax
            x0=x-xmin
            y0=y-ymin
            I2=s(x,y)
            if(I2>mint)THEN
               k=k+1
               Sum=sum+I2
            endif
            g(x0,y0)=s(x,y)
         end do
      end do
      
      do y0=0, dim2
         do x0=0, dim1
            WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
            !23        FORMAT(F10.4,AD VANCE='NO')      if you try to use a format line, 
            !                                           you're gonna have a bad time...
         end do
         Write(j,*)
      end do
      
      
      sum=sum/k
      WRITE(z,'(F10.4,i3)')Sum,i
      WRITE(j,*)
      
      CLOSE(unit=i)
      close(j)
   enddo loop1
   
   print *,'end of do loop 1 of 3'
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   ! SECOND ITERATION OF THE LOOP             !
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   
   
   
   loop2: do i = 10 ,99 
      j=i+400
      sum =0
      WRITE(cnum2,'(i2)')i
      n2=n1//n6//cnum2//n3
      n4=n10//posx//posy//u1//n6//cnum2//n3  !41+6+4+3=55
      
      ! print*, n4
      open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
      Open(unit=j,file=n4,STATUS='UNKNOWN',ACTION='WRITE')
      READ(i,*)s
      
      k=0
      do y=ymin, ymax
         do x=xmin, xmax
            x0=x-xmin
            y0=y-ymin
            I2=s(x,y)
            if(I2>mint)THEN
               k=k+1
               Sum=sum+I2
            endif
            g(x0,y0)=s(x,y)
         end do
      end do
      
      
      do y0=0,dim2
         do x0=0,dim1          
            WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
            !super cereal here (see loop 3)
         end do
         WRITE(j,*)
      end do
      WRITE(j,*)
      
      !I've never liked the print command
      !to print is human, to write divine
      !write(z,'("This shall be simple:",f10.4)',ADVANCE='NO')demand
      !I wonder if anybody is ever going to read this...
      
      sum=sum/k
      WRITE(z,'(F10.4,I3)')sum,i 
      WRITE(j,*)
      
      CLOSE(unit=i)
      CLOSE(j)
   enddo loop2
   print  *,"ending the second loop"
   !----------------------------------!
   ! THIRD ITERATION OF THE LOOP      !
   !----------------------------------!
   
   
   loop3: do i = 100 ,360 !in one frame
      j=i+400
      sum =0
      WRITE(cnum3,'(i3)') i
      n2=n1//cnum3//n3
      n4=n10//posx//posy//u1//cnum3//n3  !41+6+4+3=55
      !   print*, n4
      open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
      Open(unit=j,file=n4,STATUS='UNKNOWN',ACTION='WRITE')
      READ(i,*)s
      K=0
      do y=ymin, ymax
         do x=xmin, xmax
            x0=x-xmin
            y0=y-ymin
            I2=s(x,y)
            if(I2>mint)THEN
               Sum=sum+I2
               K=K+1
            endif
            g(x0,y0)=s(x,y)
         end do
      end do
      
      do y0=0,dim2
         do x0=0,dim1
            WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
            !I'm super cereal when I say this needs to be done in one line.
            ! I'm also super cerial when I say I don't know why....
         end do
         WRITE(j,*)
      end do
      WRITE(j,*)
      
      sum=sum/K
      WRITE(z,'(F10.4," ",i3)')sum,i
      !44  FORMAT(F10.4,i3)
      CLOSE(unit=i)
      CLOSE(j)
   enddo loop3
   print *,"ending the 3rd loop"
   
   CLOSE(23) 
   !in a nutshell every time we found a peak this ran, we opened all the files, read the values near the peak in a box the size of xmax-xmin by ymax-ymin and wrote these values along with their intensity to a file somewhere....
END DO

Print *,'You are now exiting the subroutine, Goodbye!' !
END SUBROUTINE time_ev

!|--------------------------------------~---------------------------------------------------------------|
!| 				This Space							        |
!|					Intentionally						        |
!|				Left Blank							        |
!|------------------------------------------------------------------------------------------------------|


SUBROUTINE Gauss(TABLER,TABLEI)

IMPLICIT NONE
INTEGER,DIMENSION(0:6,0:60) :: TABLEI
REAL,DIMENSION(0:4,0:60) :: TABLER
INTEGER :: xmin,ymin,xmax,ymax,i,xpeak,ypeak,rup,rdown,r,j,spots,x,y,rmax,k,zztop,abigo
REAL :: r0,Intensity,r1,r2,x0,y0,x1,y1,PI,circ
REAL, DIMENSION(0:299,0:299) :: s,q,p

! Used to put variables in file names
CHARACTER(len=1) :: cnum ! 1-9
CHARACTER(len=2) :: cnum2 ! 10-99
CHARACTER(len=3) :: cnum3 !100-360
CHARACTER(len=4) :: posx !xvalue initially found at  --\___might want to change
CHARACTER(len=4) :: posy !y value initially found at --/    to peak center

! Used to mark directories

CHARACTER (len=33):: N1 !/HOME/MICHAEL/RESEARCH/DJ11/DJ11F+ OTHER
CHARACTER(len=34) :: nu 
CHARACTER (len=38):: n2 !output path same as input+ other
CHARACTER (len=39):: n2a !output path same as input+ other
CHARACTER (len=40):: n2b !output path same as input+ other
CHARACTER(len=4) :: n3 !.dat
CHARACTER(len=1) :: n6 !used to add 0
CHARACTER(len=1) :: n7 ! used to add 0
CHARACTER(len=1) :: u1 ! used to add _
CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=48):: n5 !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=49):: n5a !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=50):: n5b !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 


!this will ultimately turn the 2d spots into 1d over the radial distance. Basically we average all phi and wind up with a 1d Gaussian.
OPEN(unit=901,FILE='1DGauss.data',STATUS='UNKNOWN',ACTION='WRITE') ! This will be the output file. 
PRINT *,"Please indicate the maximum desired radius from the center of the peak"
READ (*,*)rmax

n3='.dat'
u1='_'
n1='/home/ramalab/research/dj11/dj11f'
nu='/home/ramalab/research/Gauss/dj11G' !The G is for Gauss, or Gangster: You decide
PI=3.14159

Do abigo=0,60
   xmin=TABLEI(3,abigo)
   ymin=TABLEI(4,abigo)
   xmax=TABLEI(5,abigo)
   ymax=TABLEI(6,abigo)
   xpeak=TABLEI(1,abigo) 
   ypeak=TABLEI(2,abigo)

   write(posx,'(i4)')xpeak
   write(posy,'(i4)')ypeak

   ! LOOP 1 OF 3 !

   do i=1,9
      k=i+400
      WRITE(CNUM,'(i1)')i

      n2=n1//cnum//n3 !33+5 =38 
      n5=nu//posx//posy//u1//cnum//n3 !34+7+1+4
      
      OPEN(unit=k,file=n5,STATUS='NEW',ACTION='WRITE')
      open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
      READ(i,*)s
      
      !defining out coordinate system
      x0=xpeak
      y0=ypeak
      r0=sqrt(x0*x0+y0*y0)
      rup=1
      rdown=0
      
      do r=0,rmax
         zztop=0
         Intensity=0.0
         rup=r+1
         do x=xmin,xmax
            do y=ymin,ymax
               x1=x
               y1=y
               r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))  
               if(r2<rup.AND.r2>=r)THEN
                  Intensity=Intensity+s(x,y)
                  zztop=zztop+1
               endif
               
            end do
         end do
         Intensity=Intensity/zztop
 !        if(r>0)THEN
 !           Circ=2*PI*r
 !           Intensity=Intensity/Circ
 !        endif
         WRITE(k,'(F12.4)')Intensity
      end do
      CLOSE(i)
      CLOSE(k)
   end do
   
!LOOP 2 OF 3
   do i=10,99

      k=i+400
      WRITE(CNUM2,'(i2)')i

      n2a=n1//cnum2//n3 !33+7 =40 
      n5a=nu//posx//posy//u1//cnum2//n3 !34+7 =41 

      OPEN(unit=k,file=n5a,STATUS='NEW',ACTION='WRITE')
      open(unit=i,file=n2a,STATUS='OLD',ACTION='READ')
      READ(i,*)s

      !defining out coordinate system
      x0=xpeak
      y0=ypeak
      r0=sqrt(x0*x0+y0*y0)
      rup=1
      rdown=0
      
      do r=0,rmax
         zztop=0
         rup=r+1
         Intensity=0.0
         do x=xmin,xmax
            do y=ymin,ymax
               x1=x
               y1=y
               r1=sqrt(x1*x1+y1*y1)
               r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))       

               if(r2<rup.AND.r2>=r)THEN
                  Intensity=Intensity+s(x,y)
                  zztop=zztop+1
               endif
               
            end do
         end do
         Intensity=Intensity/zztop
!         if(r>0)THEN
!            Circ=2*pi*r
!            Intensity=(Intensity/Circ)
!         endif
         WRITE(k,'(F12.4)')Intensity
      end do
      CLOSE(i)
      CLOSE(k)
   end do  
   
! LOOP 3 OF 3
   do i=100,360

      k=i+400
      WRITE(CNUM3,'(i3)')i

      n2b=n1//cnum3//n3 !33+7 =40 
      n5b=nu//posx//posy//u1//cnum3//n3 !34+7 =41 

      open(unit=i,file=n2b,STATUS='OLD',ACTION='READ')
      OPEN(unit=k,file=n5b,STATUS='NEW',ACTION='WRITE')

      READ(i,*)s
      !defining out coordinate system

      x0=xpeak
      y0=ypeak
      r0=sqrt(x0*x0+y0*y0)
      rup=1
      rdown=0

      do r=0,rmax
         j=0
         zztop=0
         Intensity=0.0
         rup=r+1
         do x=xmin,xmax
            do y=ymin,ymax
               x1=x
               y1=y
               r1=sqrt(x1*x1+y1*y1)
               r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))  

               if(r2<rup.AND.r2>=r)THEN
                  Intensity=Intensity+s(x,y)
                  zztop=zztop+1
               endif
               
            end do
         end do
         Intensity=Intensity/zztop
!         if(r>0)THEN
!            Circ=2*pi*r
!            Intensity=(Intensity/Circ)
!         endif
         WRITE(k,'(F12.4)')Intensity
      end do
      close(i)
      close(k)
   end do
END DO

END SUBROUTINE Gauss


!SUBROUTINE PARAMETER(mint,rawin,tevout,home)
!
!INTEGER ::val,vala,valb
!
!235 Print *,"|---------------------------------------------------------------------|"
!Print *,"|      welcome, please select an option from the list below           |"
!Print *,"|---------------------------------------------------------------------|"
!Print *,"|      0 -- Program Constants                                         |"
!print *,"|      1 -- File Directories                                          |"
!print *,"|      2 --                                                           |"
!print *,"|      2 --                                                           |"
!print *,"|      2 --                                                           |"
!print *,"|      5 --                                                           |"
!print *,"|      6 -- Exit to main menu                                         |"
!print *,"|_____________________________________________________________________|"
!Read (*,*) val


!if(val.eq.0)THEN
!234 Print *,"|---------------------------------------------------------------------|"
!   Print *,"|      What Constant would you like to fiddle with?                   |"
!   Print *,"|---------------------------------------------------------------------|"
!   Print *,"|      0 -- 'MINT' the minimum threshold                              |"
!   print *,"|      1 -- 'ARRAY SIZE' The size of the imeage input                 |"
!   print *,"|      2 --  NA                                                       |"
!   print *,"|      2 --  NA                                                       |"
!   print *,"|      2 --  NA                                                       |"
!   print *,"|      5 --  NA                                                       |"
!   print *,"|      6 -- Exit to main menu                                         |"
!   print *,"|_____________________________________________________________________|"
!   Read (*,*) vala
!   IF(vala.eq.0)THEN
!      Print*,"The current Minimum threshold is: ",mint
!      print*,"Please enter the desired minimum threshold"
!      READ(*,*)mint
!      print*,"your value has been changed"
!      goto 243
!   elseif(valva.eq.1)THEN 
!      print*,"The current array size is: ",ARRAY_SIZE
!      print*,"Please indicate the desired image size: "
!      print*,'NOTE: The array begins at zero and therefore array size is the number of pixels minus one.'
!      READ(*,*)ARRAY_SIZE
!   elseif(vala.eq.6)
!   endif
   
!else if(val.eq.1) THEN
!233 Print *,"|---------------------------------------------------------------------|"
!Print *,"|      Directories are tricky this menu is a long time in the making. |"
!Print *,"|---------------------------------------------------------------------|"
!Print *,"|      0 -- Raw Data Input Files                                      |"
!print *,"|      1 -- Time Evolution Output                                     |"
!print *,"|      2 -- Gaussian parameter                                        |"
!print *,"|      3 -- Home Directory (The highest common directory)             |"
!print *,"|      4 --                                                           |"
!print *,"|      5 --                                                           |"
!print *,"|      6 -- Exit to main menu                                         |"
!print *,"|_____________________________________________________________________|"
!Read (*,*) valb
!if(valb.eq.0)THEN
!elseif(valb.eq.1)THEN
!   PRINT*,"Not sure what to say here: ",rawin
!elseif(valb.eq.2)THEN
!   PRINT*,"funny you think this text is easy to write: ",tevout
!elseif(valb.eq.3)THEN
!   PRINT*,"Current home directory is: "home
!elseif(valb.eq.6)THEN
!endif
!elseif(val.eq.6)THEN
!endif


!That should send us back to main by default
!
!
!
!
!
!
!
!END SUBROUTINE PARAMETER


SUBROUTINE GAUSS_SUBTRACT(TABLER,TABLEI)
!SAME OLD VARIABLES, NEW SONG!
integer :: i,j,l,m,z,nums,spots,t,mjd,ljd,anyvariable,max,qq
CHARACTER(len=3) :: cnum
CHARACTER(len=2) :: cnum2
CHARACTER(len=1) :: cnum1
CHARACTER(len=3) :: posx
CHARACTER(len=3) :: posy
CHARACTER(len=4) :: n3
CHARACTER(len=33) :: n1,daftpunk
CHARACTER(len=1) :: n6
CHARACTER(len=1) :: n7
CHARACTER(len=1) :: u1
CHARACTER(len=44) :: mp
CHARACTER(len=40) :: nw

REAL, DIMENSION(0:299,0:299) :: S,G,G2,F
INTEGER :: XMAX,YMAX,XMIN,YMIN,XII,YII,x0,y0,dim1,dim2,k,x,y,xpeak,ypeak,responsum
REAL ::  Sum,I2,mint, fxy,x1,y1,xpeak1,ypeak1,sigma,mew
LOGICAL :: hortari

print*,"Please enter number of spots you wish to regress."
Read(*,*)max
Print*,"Would you like your values input from HEADER.dat? Note: 1 for yes and 2 for no will suffice."
Read*,responsum
if (responsum.eq.1)THEN
   OPEN(623,FILE='HEADER.dat') ! LETS be versitile here let's construct this file in the following way
   hortari=.TRUE.
elseif(responsum.eq.0)THEN
   continue
   hortari=.FALSE.
else
   Print*,"invalid response, You are now forced to manually input the values."
   hortari=.FALSE.
endif
n3='.dat'
n1='/home/ramalab/research/dj11/dj11f' !5+8+9+5+6=
daftpunk='/home/ramalab/research/regressed/'!5+8+9+11
u1='_'

! time |xpeak |ypeak | Amplitude | sigma |mew | kai /\ this way we can construct an ifloop to read one line at atime until we have the desired peak at the desired time. The file will be long but more efficent than  header for each spot, or each time.

mjd=0
ljd=0
do anyvariable=0,max
   qq=anyvariable
   !print*,"Please input t, xpeak, ypeak, amp, sigma, and mew"
   IF(hortari.eqv..TRUE.)THEN
      READ(623,FMT='(I3,1x,I3,1x,I3,1x,F6.4,1x,F6.4,1x,F6.4,1x,F6.4,1x)')t,xpeak,ypeak,Amp,Sigma,Mew
      Print *,t,xpeak,ypeak,amp,sigma,mew
   else
      PRINT*," PLEASE input the values for time, xpeak, ypeak, Amplitude, Sigma and Mew manually now: "
      READ *,t,xpeak,ypeak,Amp,Sigma,Mew
   endif
   ! REAL, DIMENSION(0:dim1,0:dim2) :: V

   WRITE(posx,'(i3)')xmin
   WRITE(posy,'(i3)')ymin
   xmin=xpeak-5
   xmax=xpeak+5
   ymin=ypeak-5
   ymax=ypeak+5
   xpeak1=xpeak
   ypeak1=ypeak

   write(posx,'(i3)')xpeak
   write(posy,'(i3)')ypeak

   if(t<10)THEN
      WRITE(n6,'(i1)')mjd
      WRITE(n7,'(i1)')ljd
      write(cnum1,'(i1)')t
      nw=n1//n6//n7//cnum1//n3 !nw=33+1+1+1+4=40    n1=33
      mp=daftpunk//posx//u1//posy//u1//n6//n7//cnum1 !33+3+3+1+1+1+1+1=5+6+33=44
   else if(t>=10.AND.t<100)THEN
      WRITE(n6,'(i1)')mjd
      WRITE(cnum2,'(i2)')t
      nw=n1//n6//cnum2//n3
      mp=daftpunk//posx//u1//posy//u1//n6//cnum2
   else
      WRITE(cnum,'(i3)')t
      nw=n1//cnum//n3
      mp=daftpunk//posx//u1//posy//u1//cnum
   endif

 !  print *,nw

   OPEN(678,FILE=mp,STATUS='UNKNOWN',ACTION='WRITE')
   OPEN(qq,FILE=nw,STATUS='OLD',ACTION='READ')
   READ(qq,*)S

   !print*,xmin,ymin
   G=0
   do x=xmin,xmax
      do y=ymin,ymax
         !print*,"IN DO"
         x1=x
         y1=y
         r2=sqrt((x1*x1)+(y1*y1)) !not a vector
         r0=sqrt((xpeak1*xpeak1)+(ypeak1*ypeak1)) !still not a vector
         r=sqrt((x1-xpeak)*(x1-xpeak)+(y1-ypeak)*(y1-ypeak))
         f(x,y)=A*EXP(((r-mew)*(r-mew))/(sigma*sigma)) ! we can fix that later UBER PROBLEM! A- THESE AREN"T VECTORS 
         G(x,y)=f(x,y)
      end do
   end do
  ! print *," OUT DO"
 
   !  S=S2 !why did I do this? 
   !Print*,"DO THE DEW" ! we get this far befor ethe segfault?
   do  x=0,299
      do y=0,299
         G2(x,y)=S(x,y)-G(x,y) ! we need original-G
      end do
   end do
   PRINT*,"Between Two loops with zack galifinakis"
   do x=0,299
      do y=0,299
         WRITE(678,'(F12.4)',ADVANCE='NO')G2(x,y)
      end do
      Write(678,*)
   end do
   !print*,"Dew did that do?"
CLOSE(qq)
CLOSE(678)

end do
!print*,"jewjitdew?"
END SUBROUTINE GAUSS_SUBTRACT



! SUBROUTINE CONVERGE







!END SUBROUTINE CONVERGE


!This is written to be a possible subroutine added so that we can set up our alpha tables






SUBROUTINE ARPHA(tabler,TABLEI)
  
  IMPLICIT NONE 

  ! we need to find and index the peaks
  !these variables are for use with the table array and should not be used in any other array/loop
  REAL :: gnum,xpeak,ypeak,phipeak,rpeak
  REAL, DIMENSION(0:4,0:60) :: tableR
  INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  INTEGER :: i,x,y,xmin,xmax,ymin,ymax,k
  
  
  !these are defined for everything else
  INTEGER p,qq,q,pq,qp,gs,g,gn,gg,ggp,gp,zz,SIZE,gngn,gr,gsgs,dim,LB,UB,BIN
  Real :: SUM, SUM2, SUM3, alphat, avgalpha,RISE,RUN,r0,x0,y0,INT,alphaerr,SIGMA,SIGMA1
  REAL,SAVE :: MINT=50.0
  REAL, DIMENSION(0:299,0:299) :: S
  REAL, DIMENSION(0:99) :: alpha
  
  CHARACTER(len=3) :: cnum
  CHARACTER(len=4) :: dat
  CHARACTER(len=33) :: n1
  CHARACTER(len=40) :: n2
  PRINT *,"Please enter the lower bound, upper bound, and number of bins: "
  READ(*,*)LB,UB,BIN

  k=0
  N1='/home/ramalab/research/dj11/dj11f'
  x0=142.0
  y0=169.0
  r0=sqrt(x0*x0+y0*y0)
  dat='.dat'
 ! MINT=80.0
  do i=0,60

     gnum=tablei(0,k)
     xpeak=tablei(1,k)
     ypeak=tablei(2,k)
     xmin=tablei(3,k)
     ymin=tablei(4,k)
     xmax=tablei(5,k)
     ymax=tablei(6,k)
 !    PRINT *,gnum,xpeak,ypeak
     gsgs=0
     gngn=BIN-1
     k=k+1
     do g=gsgs,gngn !this loop is strait forward to me before making changes I recommend sending me and email.
        gr=LB+5*g
        gg=gr+5
        ggp=UB-1*(gr-LB)
        gp=ggp-5
        SUM2=0
        SUM3=0
    
!        PRINT *,"LOOP1 here I come"
        
        Loop1: Do p=gr,gg

           WRITE(CNUM,'(i3)')p
           n2=n1//cnum//dat !33+7 =40 
           open(unit=23,file=n2,STATUS='OLD',ACTION='READ')
           READ(23,*)S
           dim=0
           SUM=0
           do y=ymin,ymax		
              do x=xmin,xmax
                 INT=S(X,Y)
                 IF(INT>MINT)THEN
                    SUM=SUM+S(x,y)
                    DIM=DIM+1
                 ENDIF
              end do
           end do
           SUM=SUM/DIM	
           Sum2=SUM2+SUM     
           close(23)
        ENDDO LOOP1
        
 !       PRINT *,"Loop 2 looking at you"
        Loop2: Do p=gp,ggp
           SUM=0
           WRITE(CNUM,'(i3)')p
           n2=n1//cnum//dat !33+7 =40 
           open(unit=43,file=n2,STATUS='OLD',ACTION='READ')
           READ(43,*)S
           dim=0
           do y=ymin,ymax		
              do x=xmin,xmax
                 INT=S(X,Y)
                 IF(INT>MINT)THEN
  !                  print*,S(x,y)
                    SUM=SUM+S(x,y)	
                    DIM=DIM+1
                 ENDIF
             end do
           end do
 
           SUM=SUM/DIM
           Sum3=Sum3+sum	
           close(43)
        ENDDO LOOP2
  !      PRINT *,"ZZ what you got for me"
        zz=g
        Rise=Sum3-Sum2
        Run=ggp-gg
        alpha(zz)=Rise/(5*Run)

  !      PRINT *,alpha(zz)
         
     enddo
     !    PRINT *,"end do now to"
     alphat=0
     do zz=0,BIN-1
        alphat=ALPHAT+alpha(zz)
     end do
     avgalpha=log(alphat)/log(BIN)
 !    print *,avgalpha
     TABLER(2,i)=avgalpha
     ! Error Calculation
     SIGMA1=0
     do zz=0,BIN-1
        SIGMA1=(log(alpha(zz))-log(avgalpha))*(log(alpha(zz))-log(avgalpha))+SIGMA1
     end do
     SIGMA=SIGMA1/BIN
     alphaerr=sqrt(SIGMA)
     TABLER(4,i)=alphaerr
  end do
END SUBROUTINE ARPHA








SUBROUTINE WRIT (TABLER,TABLEI)
REAL, DIMENSION(0:4,0:60) :: TABLER
INTEGER, DIMENSION(0:6,0:60) :: TABLEI
INTEGER x,y
OPEN(23,FILE='TABLE.dat',STATUS='NEW',ACTION='WRITE')

do y=0,60
   do x=0,6
      WRITE(23,'(I4)', ADVANCE='NO')tableI(x,y)
   end do
   do x=0,4
      WRITE(23,'(F10.4)',ADVANCE='NO')tableR(x,y)
   end do
   WRITE(23,*)
end do
CLOSE(23)
END SUBROUTINE WRIT




SUBROUTINE IMPORT(TABLER,TABLEI)

REAL, DIMENSION(0:4,0:60) :: TABLER
INTEGER, DIMENSION(0:6,0:60) :: TABLEI
INTEGER x,y,y5

OPEN(23,FILE='TABLE.dat',STATUS='OLD',ACTION='READ')

do y=0,60
   READ(23,'(7I4)',ADVANCE='NO')tableI(0,y),tableI(1,y),tableI(2,y),tableI(3,y),tableI(4,y),tableI(5,y),tableI(6,y)
   READ(23,'(5F10.4)')tableR(0,y),tableR(1,y),tableR(2,y),tableR(3,y),tableR(4,y)
end do
CLOSE(23)
END SUBROUTINE IMPORT

SUBROUTINE SORT(TABLER,TABLEI)

REAL, DIMENSION(0:4,0:60) :: TABLER, TEMP
INTEGER, DIMENSION(0:6,0:60) :: TABLEI
INTEGER x,y

Real :: Phi0, dphi, Rpeak, Phipeak,r0,x1,y1,r1,r2,ref,radangle,minang,maxang,xmin,ymin,dphi1 
Integer :: i, xpoint, ypoint, grp,y5, xmi, ymi
REAL,SAVE :: x0=142,y0=169

PRINT*,"PLEASE input the angle for spots to be sorted by, in radians."
READ(*,*)radangle

r0=sqrt(x0*x0+y0*y0)

do y=0,60
   xpoint=TABLEI(1,y)
   ypoint=TABLEI(2,y)
   x1=xpoint
   y1=ypoint
   xmi=TABLEI(3,y)
   ymi=TABLEI(6,y)
   xmin=xmi
   ymax=ymi
   r1=sqrt(x1*x1+y1*y1)
   r2=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
   TABLER(1,y)=r2
   phi=acos(abs(abs(x1-x0)/r2))
   TABLER(0,y)=phi
   dphi=abs(phi-acos(abs(x1-x0)/(sqrt((xmin-x0)*(xmin-x0)+(ymax-y0)*(ymax-y0)))))
   TABLER(3,y)=dphi
end do
do y=0,40
   phi0=TABLER(0,y)
   grp=TABLEI(0,y)
   dphi=TABLER(3,y)
   do y5=0,40
      phi1=TABLER(0,y5)
      dphi1=TABLER(3,y5)
      delta=abs(phi0-phi1)
      maxang=radangle+abs(dphi/2)+abs(dphi1/2)
      minang=radangle-abs(dphi/2)-abs(dphi1/2)
      if(delta<maxang.AND.delta>minang)THEN
         TABLEI(0,y5)=grp
      else
      endif
   enddo
enddo
END SUBROUTINE SORT


SUBROUTINE Q_SPACE


IMPLICIT NONE
INTEGER,DIMENSION(0:6,0:60) :: TABLEI
REAL,DIMENSION(0:4,0:60) :: TABLER
INTEGER :: xmin,ymin,xmax,ymax,i,xpeak,ypeak,rup,rdown,r,j,spots,x,y,rmax,nums,k,zztop,circumd,circ,circum
REAL :: r0,Intensity,r1,r2,x0,y0,x1,y1,PI,dphi,phiup,phimax,phimin,phi,rmin
REAL, DIMENSION(0:299,0:299) :: s,q,p

! Used to put variables in file names
CHARACTER(len=1) :: cnum ! 1-9
CHARACTER(len=2) :: cnum2 ! 10-99
CHARACTER(len=3) :: cnum3,PHIW !100-360
CHARACTER(len=4) :: posx !xvalue initially found at  --\___might want to change
CHARACTER(len=4) :: posy !y value initially found at --/    to peak center

! Used to mark directories

CHARACTER (len=33):: N1 !/HOME/MICHAEL/RESEARCH/DJ11/DJ11F+ OTHER
CHARACTER(len=34) :: nu 
CHARACTER (len=38):: n2 !output path same as input+ other
CHARACTER (len=39):: n2a !output path same as input+ other
CHARACTER (len=40):: n2b !output path same as input+ other
CHARACTER(len=4) :: n3 !.dat
CHARACTER(len=1) :: n6 !used to add 0
CHARACTER(len=1) :: n7 ! used to add 0
CHARACTER(len=1) :: u1 ! used to add _
CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=52):: n5 !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=51):: n5a !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 
CHARACTER (len=54):: n5b !Complete OUTPUT path n1+ cnum and n3 +1 for _ added 


!this will ultimately turn the 2d spots into 1d over the radial distance. Basically we average all phi and wind up with a 1d Gaussian.
OPEN(unit=900,FILE='OUTPUT.data',STATUS='OLD',ACTION='READ')
OPEN(unit=901,FILE='1DGauss.data',STATUS='UNKNOWN',ACTION='WRITE') ! This will be the output file. 
Print *,"Please indicate the number of spots in the file. NOTE: there is no direct input functionality at this time."
Read (*,*) nums
PRINT *,"Please indicate the maximum desired radius from the center of the peak"
READ (*,*)rmax
PRINT*,"Please indicate the number of phi bins to be used"
READ(*,*)circumd !from the latin for circle, circumdare

n3='.dat'
u1='_'
n1='/home/ramalab/research/dj11/dj11f'
nu='/home/ramalab/research/Q_spc/dj11Q' 
PI=3.14159
x0=142
y0=169


Do spots=1,nums
   Read (900,'(6I4)')xmin,ymin,xmax,ymax, xpeak, ypeak
   
   write(posx,'(i4)')xpeak
   write(posy,'(i4)')ypeak

   dphi=abs(acos(sqrt((xmax-x0)*(xmax-x0)+(ymin-y0)*(ymin-y0))/x0)-acos(sqrt((xmin-x0)*(xmin-x0)+(ymax-y0)*(ymax-y0))/x0))/circumd
   phimin=acos(sqrt((xmax-x0)*(xmax-x0)+(ymin-y0)*(ymin-y0))/x0)
   phimax=acos(sqrt((xmin-x0)*(xmin-x0)+(ymax-y0)*(ymax-y0))/x0)
   R0=sqrt((xpeak-x0)*(xpeak-x0)+(ypeak-y0)*(ypeak-y0))
   ! LOOP 1 OF 3 !
   circum=circumd-1



      do i=1,9
         !defining out coordinate system
         do circ=0,circum
            phiup=phimin+dphi*(circ+1)
            k=i+400
            WRITE(CNUM,'(i1)')i
            WRITE(PHIW,'(F3.0)')phiup
            n2=n1//cnum//n3 !33+5 =38 
            n5=nu//posx//posy//u1//cnum//u1//PHIW//n3 !34+7+1+4
            
            OPEN(unit=k,file=n5,STATUS='new',ACTION='WRITE')
            open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
            READ(i,*)s
            do r=0,rmax  
               rmin=r
               do x=xmin,xmax
                  do y=ymin,ymax
                     x1=x
                     y1=y
                     phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
                     r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
                     if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
                        Intensity=S(x,y)+Intensity
                     endif
                  end do
               end do
            
               Intensity=Intensity/2
               !        if(r>0)THEN
               !           Circ=2*PI*r
               !           Intensity=Intensity/Circ
               !        endif 
               WRITE(k,'(F12.4," ",I3)')Intensity,r
            end do
            CLOSE(i)
            CLOSE(k)
         enddo
      enddo
      !LOOP 2 OF 3
      do i=10,99
         do circ=0,circum
            phiup=phimin+dphi*(circ+1)
            k=i+400
            WRITE(CNUM2,'(i2)')i
            WRITE(PHIW,'(F3.0)')phiup
            n2a=n1//cnum2//n3 !33+7 =40 
            n5a=nu//posx//posy//u1//cnum2//n3 !34+7 =44 
            
            OPEN(unit=k,file=n5a,STATUS='NEW',ACTION='WRITE')
            open(unit=i,file=n2a,STATUS='OLD',ACTION='READ')
            READ(i,*)s

            do r=0,rmax  
               rmin=r
               do x=xmin,xmax
                  do y=ymin,ymax
                     x1=x
                     y1=y
                     phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
                     r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
                     if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
                        Intensity=S(x,y)+Intensity
                     endif
                  end do
               end do
            
               Intensity=Intensity/2
               !        if(r>0)THEN
               !           Circ=2*PI*r
               !           Intensity=Intensity/Circ
               !        endif 
               WRITE(k,'(F12.4," ",I3)')Intensity,r
            end do
            
            CLOSE(i)
            CLOSE(k)
         enddo
      enddo
      
      ! LOOP 3 OF 3
      do i=100,360
         
         !defining out coordinate system
         do circ=0,circum
            phiup=phimin+dphi*(circ+1)
            k=i+400
            WRITE(CNUM3,'(i3)')i
            WRITE(PHIW,'(F3.0)')phiup
            n2b=n1//cnum3//n3 !33+7 =40 
            n5b=nu//posx//posy//u1//cnum3//u1//phiw//n3 !34+7 =44 
            
            open(unit=i,file=n2b,STATUS='OLD',ACTION='READ')
            OPEN(unit=k,file=n5b,STATUS='NEW',ACTION='WRITE')
            do r=0,rmax  
               rmin=r
               do x=xmin,xmax
                  do y=ymin,ymax
                     x1=x
                     y1=y
                     phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
                     r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
                     if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
                        Intensity=S(x,y)+Intensity
                     endif
                  end do
               end do
               Intensity=Intensity/2
               !        if(r>0)THEN
               !           Circ=2*PI*r
               !           Intensity=Intensity/Circ
               !        endif 
               WRITE(k,'(F12.4," ",I3)')Intensity,r
            end do
            CLOSE(i)
            CLOSE(k)
         enddo
      ENDDO
   enddo

END SUBROUTINE Q_SPACE