!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 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.> !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(Ixmax)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 (IImax) THEN xpeak=x0 ypeak=y0 end if if(I2ymax)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(I22.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=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=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=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(deltaminang)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(phiphimin.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(phiphimin.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(phiphimin.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