peakmod.f90

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! 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