peakmod.f90


SUBMITTED BY: Guest

DATE: April 11, 2013, 1:14 a.m.

FORMAT: Text only

SIZE: 48.9 kB

HITS: 838

  1. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2. ! Peak.f90 v1.0.1 !
  3. ! Summer 2012 !
  4. ! note: functionality has been added several times !
  5. ! !
  6. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  7. ! Release Notes: ! !
  8. ! -fixed issues with minimum threshold. Anything below about 100 !
  9. ! is unsuitble !
  10. ! " Gauss subtract: !
  11. ! Says end of file in RTE no idea why !
  12. ! " Parameters !
  13. ! -Very buggy but not important and bottom of the priority list !
  14. ! " measure peak now measures alpha value !
  15. ! " write... writes tons of stuff to a file, tis good !
  16. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  17. PROGRAM MENU
  18. Integer :: val
  19. LOGICAL :: C
  20. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  21. REAL, DIMENSION(0:4,0:60) :: TABLER
  22. C=.TRUE.
  23. Tablei=0
  24. tabler=0
  25. DO WHILE (C.EQV..TRUE.)
  26. Print *,"|---------------------------------------------------------------------|"
  27. Print *,"| welcome to SAT v 1.0 |"
  28. print *,"| please select an option from the list below |"
  29. Print *,"|---------------------------------------------------------------------|"
  30. Print *,"| 0 -- Find new peaks |"
  31. print *,"| 1 -- Measure Peaks growth |"
  32. print *,"| 2 -- Prepare 1D Gaussian Bins |"
  33. print *,"| 3 -- Prepare peaks for imaging |"
  34. print *,"| 4 -- Subtract Gaussian from Data [run-time-FDNE] |"
  35. print *,"| 5 -- Write Table values to File |"
  36. print *,"| 6 -- Import Table values from File |"
  37. print *,"| 7 -- SORT Alpha Table |"
  38. print *,"| 8 -- Q Space Gaussian-- pretty buggy |"
  39. print *,"| 9 -- PARAMETERS N/A |"
  40. print *,"| 10 -- exit |"
  41. print *,"|_____________________________________________________________________|"
  42. Read (*,*) val
  43. If (val.eq.0) THEN
  44. call Peakfinder(TABLER,TABLEI)
  45. !goto 23
  46. else If (val.eq.1) THEN
  47. Call ARPHA(tabler,tablei)
  48. ! goto 23
  49. else If (val.eq.2) THEN
  50. call Gauss(TABLER,TABLEI)
  51. else If (val.eq.3) THEN
  52. call time_ev(TABLER,TABLEI)
  53. ! goto 23!
  54. else If (val.eq.4) THEN
  55. call Gauss_subtract(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
  56. ! goto 23
  57. else IF(val.eq.5)THEN
  58. call WRIT(TABLER,TABLEI)
  59. else If (val.eq.6) THEN
  60. call IMPORT(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
  61. else If (val.eq.7) THEN
  62. call SORT(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
  63. ! goto 23
  64. else If (val.eq.8) THEN
  65. call Q_SPACE(TABLER,TABLEI) ! oncve we have the main functions done this can be added, perhaps even to 3
  66. ! goto 23
  67. else If (val.eq.10) THEN
  68. C=.FALSE.
  69. EXIT
  70. else
  71. print *, "please select a valid menu option"
  72. ! goto 23
  73. endif
  74. !Kenny Peou is great,
  75. !Menu would not have happened,
  76. !without his input...
  77. END DO
  78. STOP
  79. END PROGRAM MENU
  80. !Table format
  81. !TABLEI
  82. ! Group# | xpeak | ypeak | xmin |ymin |xmax |ymax|
  83. !
  84. !TABLER
  85. !phi peak \ rpeak \ alpa \ dphi
  86. !
  87. !
  88. SUBROUTINE PEAKFINDER(tabler,tablei)
  89. !IMPLICIT NONE
  90. REAL, DIMENSION(0:4,0:60) :: tabler
  91. INTEGER, DIMENSION(0:6,0:60) :: tablei
  92. REAL, DIMENSION(0:299,0:299) :: s,q,p
  93. INTEGER :: x, y, x1, y1, x2, y2, xi, yi, xf, yf, xii, yii, xmax0, x0, y0, xmax1, ymax1, xmin1, ymin1, wxpeak, wypeak
  94. LOGICAL Spot, Spot2, Newspots, Write, Gauss
  95. INTEGER,SAVE :: xmin,ymin,xmax,ymax, lambda
  96. REAL :: xpeak, ypeak, I, I2, Imax, rmax, rmin, r, mint, r0
  97. !just opening the files and assigning them to an array
  98. OPEN(923,FILE='/home/ramalab/research/dj11/dj11f355.dat',STATUS='OLD',ACTION='READ')
  99. OPEN(900,FILE='OUTPUT.data')
  100. READ(923,*)s
  101. lambda=0
  102. !defining out coordinate system
  103. !r0=sqrt(x0*x0+y0*y0)
  104. x0=142
  105. y0=169
  106. rmin=r0+82
  107. rmax=rmin+48
  108. mint=100.0
  109. Spot=.FALSE.
  110. Spot2=.FALSE.
  111. ymin=0
  112. ymax=0
  113. xmin=0
  114. xmax=0
  115. !Note this next line may need ot be moved, I am not sure that this is entirely kosher
  116. q=s
  117. p=s
  118. ! WHERE(s>rmin.and.s<rmax)
  119. !easiest way to do this loop hence why we do this in Fortran not C >.>
  120. !when we find a point in the two rings that meets this condition of a min value
  121. do x=0,299
  122. do y=0,299
  123. q(x,y)=p(x,y)
  124. I=q(x,y)
  125. if (I >= mint)then
  126. ! PRINT *,"Found one!"
  127. ! Print *,lambda
  128. tablei(0,lambda)=lambda
  129. ! lambda=lambda+1
  130. ! Print*,I,x,y
  131. Spot=.TRUE.
  132. spot2=.TRUE.
  133. !variables come
  134. !also variables go
  135. !the below are saved
  136. x1=x
  137. y1=y
  138. ymin=y
  139. xmin=x
  140. ymax=y
  141. xmax=x
  142. p=q
  143. Imax=0.0
  144. !PRint *," ENTERING ACTUAL FIRST LOOP"
  145. DO WHILE(Spot.eqv..TRUE.)
  146. I=p(x0,y0)
  147. If(I < mint)THEN
  148. Spot=.FALSE.
  149. endif
  150. y0=y0+1
  151. y0=ymax
  152. END DO
  153. !PRINT*,"ENDING FIRST LOOP"
  154. ymin1=ymin
  155. ymax1=ymax
  156. !print*,"starting loop1"
  157. Do y0=ymin1,ymax1
  158. SPOT2=.TRUE.
  159. x0=x1
  160. DO WHILE(spot2.eqv..TRUE.)
  161. I=p(x0,y0)
  162. if(I>Imax)THEN
  163. !some Haikus composed for the issue
  164. !---Here is the problem
  165. !What it is I do not see
  166. !tab is off in loop?---
  167. !---Line one-thirty five
  168. !One thirty six,between them,
  169. !Something goes awry.
  170. !------------------------new style-----------Actual Rhyme ---- Ideal
  171. !when implicit equals some, A A
  172. !values undefined, B B
  173. !are no problem of any kind. B A
  174. !though if we set implicit none, C B
  175. !then we find, B C
  176. !programming becomes more fun. C B
  177. !implicit defining creates more problems, D C
  178. !than it illuminates answers, E A
  179. !yeah, I am out of ideas, schmoblems... D C
  180. xpeak=x0
  181. ypeak=y0
  182. Imax=I
  183. endif
  184. if(I<mint)THEN
  185. Spot2=.FALSE.
  186. endif
  187. x0=x0+1
  188. if(x0>xmax)THEN
  189. xmax=x0
  190. endif
  191. end do
  192. SPOT2=.TRUE.
  193. x0=x1
  194. DO WHILE(spot2.eqv..TRUE.)
  195. I=p(x0,y0)
  196. if(I>Imax) THEN
  197. xpeak=x0
  198. ypeak=y0
  199. end if
  200. if (I<mint)THEN
  201. Spot2=.FALSE.
  202. endif
  203. x0=x0-1
  204. if(x0<xmin)THEN
  205. xmin=x0
  206. endif
  207. end do
  208. end do
  209. xmin1=xmin
  210. xmax1=xmax
  211. Do x0=xmin1,xmax1
  212. Spot2=.TRUE.
  213. y0=y1
  214. DO WHILE(spot2.eqv..TRUE.)
  215. I2=p(x0,y0)
  216. if(I2>Imax) THEN
  217. xpeak=x0
  218. ypeak=y0
  219. end if
  220. if(I2<mint)THEN
  221. Spot2=.FALSE.
  222. end if
  223. y0=y0+1
  224. if(y0>ymax)THEN
  225. ymax=y0
  226. endif
  227. end do
  228. SPOT2=.TRUE.
  229. y0=y1
  230. !Print *,"alg in process pease stand by"
  231. DO WHILE(spot2.eqv..TRUE.)
  232. I2=p(x0,y0)
  233. if(I2>Imax) THEN
  234. xpeak=x0
  235. ypeak=y0
  236. Imax=I
  237. end if
  238. IF(I2<mint)THEN
  239. Spot2=.FALSE.
  240. endif
  241. y0=y0-1
  242. if(y0<ymin)THEN
  243. ymin=y0
  244. endif
  245. !print *,ymin
  246. end do
  247. end do
  248. do x0=xmin,xmax
  249. do y0=ymin,ymax
  250. p(x0,y0)=0.0
  251. end do
  252. end do
  253. If(xmin>2.AND.ymin>2)THEN
  254. xmin=xmin-2
  255. ymin=ymin-2
  256. endif
  257. if(xmax<297.AND.ymax>297)THEN
  258. xmax=xmax+2
  259. ymax=ymax+2
  260. endif
  261. wxpeak=NINT(xpeak)
  262. wypeak=NINT(ypeak)
  263. Write(900,'(6I4)')xmin, ymin, xmax, ymax, wxpeak, wypeak
  264. !Print*,xpeak,ypeak,xmin,ymin !we get stuck like 1/3 of the way through, why?
  265. tablei(1,lambda)=xpeak
  266. tablei(2,lambda)=ypeak
  267. tablei(3,lambda)=xmin
  268. tablei(4,lambda)=ymin
  269. tablei(5,lambda)=xmax
  270. tablei(6,lambda)=ymax
  271. ! print *,tablei(0,lambda)
  272. !Late july, a Haiku on issues unknown
  273. !---While I debug this,
  274. !I begin to take notice,
  275. !How much is broken. ---
  276. !--- Array 2 9 9
  277. !Loop goes zero three hundred
  278. !Doesn't matter at all---
  279. !---If statement shouldn't end
  280. !Add else condition to end
  281. !Doesn't matter at all---
  282. !---One third the way through
  283. !The loop fails to continue,
  284. !Why on Earth one third?---
  285. !---Every time you ask?
  286. !That is the sign of a bug
  287. !does one third matter?---
  288. !---or coincidence?
  289. !hard to believe it does not.
  290. !Hard to connect dots.
  291. !~~~~~~~~~~~~~~~~~~~~~~~~~~Let's switch styles~~~~~~~~~~~~~~~~~~~~~~~~~~~
  292. !Looping over variable x,
  293. !Looping over variable y,
  294. !as long as I don't care, try exchange,
  295. !Y before repeats the same this try.
  296. ! MY STYLE NEEDS MORE WORK THAN THIS CODE
  297. q=p
  298. lambda=lambda+1
  299. elseif(I < mint)THEN
  300. !do a little dance?
  301. end if
  302. end do
  303. end do
  304. close(900)
  305. close(923)
  306. Print*,"That should be the last of them skipper!"
  307. RETURN
  308. END SUBROUTINE PEAKFINDER !(TABLER, TABLEI)
  309. SUBROUTINE time_ev(TABLER,TABLEI)
  310. IMPLICIT NONE
  311. REAL,DIMENSION(0:4,0:60) :: TABLER
  312. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  313. ! DO NOT CHANGE THESE CHANGE THE FILES TO FIT THE STRUCTURE DO NOT CHANGE THESE!
  314. integer :: i,j,l,m,z,nums,spots
  315. ! Used to put variables in file names
  316. CHARACTER(len=1) :: cnum ! 1-9
  317. CHARACTER(len=2) :: cnum2 ! 10-99
  318. CHARACTER(len=3) :: cnum3 !100-360
  319. CHARACTER(len=3) :: posx !xvalue initially found at --\___might want to change
  320. CHARACTER(len=3) :: posy !y value initially found at --/ to peak center
  321. ! Used to mark directories
  322. CHARACTER(len=33) :: n1 !/home/michael/research/dj11/dj11f+ other
  323. CHARACTER(len=41) :: n10 !/home/michael/research/dj11/Spot_ev/dj11f+ other
  324. CHARACTER(len=37) :: n11 !/home/michael/research/dj11/IvT/dj11f+ other
  325. CHARACTER(len=40) :: n2 !output path same as input+ other
  326. CHARACTER(len=4) :: n3 !.dat
  327. CHARACTER(len=1) :: n6 !used to add 0
  328. CHARACTER(len=1) :: n7 ! used to add 0
  329. CHARACTER(len=1) :: u1 ! used to add _
  330. CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added
  331. CHARACTER(len=47) :: n5 !outputs for Intensity vs time
  332. !this one is for file format -- didn't use it
  333. CHARACTER(len=3) :: daftpunk
  334. REAL, DIMENSION(0:299,0:299) :: S,g
  335. INTEGER :: x,y,XMAX,YMAX,XMIN,YMIN,XII,YII,x0,y0,dim1,dim2,k,xpeak,ypeak
  336. REAL :: Sum,I2,mint
  337. !OPENS every file, gets info for
  338. OPEN(unit=900,FILE='OUTPUT.data',STATUS='old',ACTION='read')
  339. Print *,"Please indicate the number of spots in the file."
  340. Read (*,*) nums
  341. Do spots=1,nums
  342. Read(900,'(6I3)')xmin,ymin,xmax,ymax,xpeak,ypeak
  343. ! Stuff I probably don't need
  344. X=XMIN
  345. Y=YMIN
  346. sum=0
  347. mint=100.0
  348. ! Stuff I may need?
  349. dim1=xmax-xmin
  350. dim2=ymax-ymin
  351. ! REAL, DIMENSION(0:dim1,0:dim2) :: V
  352. WRITE(posx,'(i3)')xmin
  353. WRITE(posy,'(i3)')ymin
  354. WRITE(daftpunk,'(i3)')dim1 !totally not used
  355. !Definted for all Three Loops
  356. N1='/home/ramalab/research/dj11/dj11f'
  357. N10='/home/rabalab/research/dj11/Spot_ev/dj11f'
  358. N11='/home/ramalab/research/dj11/Ivt/dj11f'
  359. n3='.dat'
  360. u1='_'
  361. n5=n11//posx//posy//n3 !37+10=47
  362. m=0
  363. l=0
  364. WRITE(n6,'(i1)')m
  365. WRITE(n7,'(i1)')l
  366. !I'm just baffled how this won't worked....
  367. OPEN(23,file=n5,STATUS='NEW',ACTION='WRITE')
  368. PRINT *,'You are now in the subroutine, enjoy it while it lasts >.>'
  369. print*,"dim1, dim2, xmin, ymin, xmax, and ymax: "
  370. print*,dim1,dim2,xmin,ymin,xmax,ymax
  371. !------------------------------------~~----------------------------------------!
  372. ! These are the loops That Try mens souls... !
  373. ! !
  374. ! First Iteration of 'THE' Loop !
  375. !------------------------------------------------------------------------------!
  376. LOOP1: DO i=1,9
  377. j=i+400
  378. !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
  379. sum=0
  380. WRITE(CNUM,'(i1)')i
  381. n2=n1//n6//n7//cnum//n3 !33+7 =40
  382. n4=n10//posx//posy//u1//n6//n7//CNUM//n3 !41+6+4+3=55
  383. !tricky way to make sure we dont have two files open with the same tag
  384. open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
  385. Open(unit=j,file=n4,STATUS='NEW',ACTION='WRITE')
  386. READ(i,*)s
  387. k=0
  388. do y=ymin, ymax
  389. do x=xmin, xmax
  390. x0=x-xmin
  391. y0=y-ymin
  392. I2=s(x,y)
  393. if(I2>mint)THEN
  394. k=k+1
  395. Sum=sum+I2
  396. endif
  397. g(x0,y0)=s(x,y)
  398. end do
  399. end do
  400. do y0=0, dim2
  401. do x0=0, dim1
  402. WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
  403. !23 FORMAT(F10.4,AD VANCE='NO') if you try to use a format line,
  404. ! you're gonna have a bad time...
  405. end do
  406. Write(j,*)
  407. end do
  408. sum=sum/k
  409. WRITE(z,'(F10.4,i3)')Sum,i
  410. WRITE(j,*)
  411. CLOSE(unit=i)
  412. close(j)
  413. enddo loop1
  414. print *,'end of do loop 1 of 3'
  415. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
  416. ! SECOND ITERATION OF THE LOOP !
  417. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
  418. loop2: do i = 10 ,99
  419. j=i+400
  420. sum =0
  421. WRITE(cnum2,'(i2)')i
  422. n2=n1//n6//cnum2//n3
  423. n4=n10//posx//posy//u1//n6//cnum2//n3 !41+6+4+3=55
  424. ! print*, n4
  425. open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
  426. Open(unit=j,file=n4,STATUS='UNKNOWN',ACTION='WRITE')
  427. READ(i,*)s
  428. k=0
  429. do y=ymin, ymax
  430. do x=xmin, xmax
  431. x0=x-xmin
  432. y0=y-ymin
  433. I2=s(x,y)
  434. if(I2>mint)THEN
  435. k=k+1
  436. Sum=sum+I2
  437. endif
  438. g(x0,y0)=s(x,y)
  439. end do
  440. end do
  441. do y0=0,dim2
  442. do x0=0,dim1
  443. WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
  444. !super cereal here (see loop 3)
  445. end do
  446. WRITE(j,*)
  447. end do
  448. WRITE(j,*)
  449. !I've never liked the print command
  450. !to print is human, to write divine
  451. !write(z,'("This shall be simple:",f10.4)',ADVANCE='NO')demand
  452. !I wonder if anybody is ever going to read this...
  453. sum=sum/k
  454. WRITE(z,'(F10.4,I3)')sum,i
  455. WRITE(j,*)
  456. CLOSE(unit=i)
  457. CLOSE(j)
  458. enddo loop2
  459. print *,"ending the second loop"
  460. !----------------------------------!
  461. ! THIRD ITERATION OF THE LOOP !
  462. !----------------------------------!
  463. loop3: do i = 100 ,360 !in one frame
  464. j=i+400
  465. sum =0
  466. WRITE(cnum3,'(i3)') i
  467. n2=n1//cnum3//n3
  468. n4=n10//posx//posy//u1//cnum3//n3 !41+6+4+3=55
  469. ! print*, n4
  470. open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
  471. Open(unit=j,file=n4,STATUS='UNKNOWN',ACTION='WRITE')
  472. READ(i,*)s
  473. K=0
  474. do y=ymin, ymax
  475. do x=xmin, xmax
  476. x0=x-xmin
  477. y0=y-ymin
  478. I2=s(x,y)
  479. if(I2>mint)THEN
  480. Sum=sum+I2
  481. K=K+1
  482. endif
  483. g(x0,y0)=s(x,y)
  484. end do
  485. end do
  486. do y0=0,dim2
  487. do x0=0,dim1
  488. WRITE(j,'(F5.1)',ADVANCE='NO')g(x0,y0)
  489. !I'm super cereal when I say this needs to be done in one line.
  490. ! I'm also super cerial when I say I don't know why....
  491. end do
  492. WRITE(j,*)
  493. end do
  494. WRITE(j,*)
  495. sum=sum/K
  496. WRITE(z,'(F10.4," ",i3)')sum,i
  497. !44 FORMAT(F10.4,i3)
  498. CLOSE(unit=i)
  499. CLOSE(j)
  500. enddo loop3
  501. print *,"ending the 3rd loop"
  502. CLOSE(23)
  503. !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....
  504. END DO
  505. Print *,'You are now exiting the subroutine, Goodbye!' !
  506. END SUBROUTINE time_ev
  507. !|--------------------------------------~---------------------------------------------------------------|
  508. !| This Space |
  509. !| Intentionally |
  510. !| Left Blank |
  511. !|------------------------------------------------------------------------------------------------------|
  512. SUBROUTINE Gauss(TABLER,TABLEI)
  513. IMPLICIT NONE
  514. INTEGER,DIMENSION(0:6,0:60) :: TABLEI
  515. REAL,DIMENSION(0:4,0:60) :: TABLER
  516. INTEGER :: xmin,ymin,xmax,ymax,i,xpeak,ypeak,rup,rdown,r,j,spots,x,y,rmax,k,zztop,abigo
  517. REAL :: r0,Intensity,r1,r2,x0,y0,x1,y1,PI,circ
  518. REAL, DIMENSION(0:299,0:299) :: s,q,p
  519. ! Used to put variables in file names
  520. CHARACTER(len=1) :: cnum ! 1-9
  521. CHARACTER(len=2) :: cnum2 ! 10-99
  522. CHARACTER(len=3) :: cnum3 !100-360
  523. CHARACTER(len=4) :: posx !xvalue initially found at --\___might want to change
  524. CHARACTER(len=4) :: posy !y value initially found at --/ to peak center
  525. ! Used to mark directories
  526. CHARACTER (len=33):: N1 !/HOME/MICHAEL/RESEARCH/DJ11/DJ11F+ OTHER
  527. CHARACTER(len=34) :: nu
  528. CHARACTER (len=38):: n2 !output path same as input+ other
  529. CHARACTER (len=39):: n2a !output path same as input+ other
  530. CHARACTER (len=40):: n2b !output path same as input+ other
  531. CHARACTER(len=4) :: n3 !.dat
  532. CHARACTER(len=1) :: n6 !used to add 0
  533. CHARACTER(len=1) :: n7 ! used to add 0
  534. CHARACTER(len=1) :: u1 ! used to add _
  535. CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added
  536. CHARACTER (len=48):: n5 !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  537. CHARACTER (len=49):: n5a !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  538. CHARACTER (len=50):: n5b !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  539. !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.
  540. OPEN(unit=901,FILE='1DGauss.data',STATUS='UNKNOWN',ACTION='WRITE') ! This will be the output file.
  541. PRINT *,"Please indicate the maximum desired radius from the center of the peak"
  542. READ (*,*)rmax
  543. n3='.dat'
  544. u1='_'
  545. n1='/home/ramalab/research/dj11/dj11f'
  546. nu='/home/ramalab/research/Gauss/dj11G' !The G is for Gauss, or Gangster: You decide
  547. PI=3.14159
  548. Do abigo=0,60
  549. xmin=TABLEI(3,abigo)
  550. ymin=TABLEI(4,abigo)
  551. xmax=TABLEI(5,abigo)
  552. ymax=TABLEI(6,abigo)
  553. xpeak=TABLEI(1,abigo)
  554. ypeak=TABLEI(2,abigo)
  555. write(posx,'(i4)')xpeak
  556. write(posy,'(i4)')ypeak
  557. ! LOOP 1 OF 3 !
  558. do i=1,9
  559. k=i+400
  560. WRITE(CNUM,'(i1)')i
  561. n2=n1//cnum//n3 !33+5 =38
  562. n5=nu//posx//posy//u1//cnum//n3 !34+7+1+4
  563. OPEN(unit=k,file=n5,STATUS='NEW',ACTION='WRITE')
  564. open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
  565. READ(i,*)s
  566. !defining out coordinate system
  567. x0=xpeak
  568. y0=ypeak
  569. r0=sqrt(x0*x0+y0*y0)
  570. rup=1
  571. rdown=0
  572. do r=0,rmax
  573. zztop=0
  574. Intensity=0.0
  575. rup=r+1
  576. do x=xmin,xmax
  577. do y=ymin,ymax
  578. x1=x
  579. y1=y
  580. r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))
  581. if(r2<rup.AND.r2>=r)THEN
  582. Intensity=Intensity+s(x,y)
  583. zztop=zztop+1
  584. endif
  585. end do
  586. end do
  587. Intensity=Intensity/zztop
  588. ! if(r>0)THEN
  589. ! Circ=2*PI*r
  590. ! Intensity=Intensity/Circ
  591. ! endif
  592. WRITE(k,'(F12.4)')Intensity
  593. end do
  594. CLOSE(i)
  595. CLOSE(k)
  596. end do
  597. !LOOP 2 OF 3
  598. do i=10,99
  599. k=i+400
  600. WRITE(CNUM2,'(i2)')i
  601. n2a=n1//cnum2//n3 !33+7 =40
  602. n5a=nu//posx//posy//u1//cnum2//n3 !34+7 =41
  603. OPEN(unit=k,file=n5a,STATUS='NEW',ACTION='WRITE')
  604. open(unit=i,file=n2a,STATUS='OLD',ACTION='READ')
  605. READ(i,*)s
  606. !defining out coordinate system
  607. x0=xpeak
  608. y0=ypeak
  609. r0=sqrt(x0*x0+y0*y0)
  610. rup=1
  611. rdown=0
  612. do r=0,rmax
  613. zztop=0
  614. rup=r+1
  615. Intensity=0.0
  616. do x=xmin,xmax
  617. do y=ymin,ymax
  618. x1=x
  619. y1=y
  620. r1=sqrt(x1*x1+y1*y1)
  621. r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))
  622. if(r2<rup.AND.r2>=r)THEN
  623. Intensity=Intensity+s(x,y)
  624. zztop=zztop+1
  625. endif
  626. end do
  627. end do
  628. Intensity=Intensity/zztop
  629. ! if(r>0)THEN
  630. ! Circ=2*pi*r
  631. ! Intensity=(Intensity/Circ)
  632. ! endif
  633. WRITE(k,'(F12.4)')Intensity
  634. end do
  635. CLOSE(i)
  636. CLOSE(k)
  637. end do
  638. ! LOOP 3 OF 3
  639. do i=100,360
  640. k=i+400
  641. WRITE(CNUM3,'(i3)')i
  642. n2b=n1//cnum3//n3 !33+7 =40
  643. n5b=nu//posx//posy//u1//cnum3//n3 !34+7 =41
  644. open(unit=i,file=n2b,STATUS='OLD',ACTION='READ')
  645. OPEN(unit=k,file=n5b,STATUS='NEW',ACTION='WRITE')
  646. READ(i,*)s
  647. !defining out coordinate system
  648. x0=xpeak
  649. y0=ypeak
  650. r0=sqrt(x0*x0+y0*y0)
  651. rup=1
  652. rdown=0
  653. do r=0,rmax
  654. j=0
  655. zztop=0
  656. Intensity=0.0
  657. rup=r+1
  658. do x=xmin,xmax
  659. do y=ymin,ymax
  660. x1=x
  661. y1=y
  662. r1=sqrt(x1*x1+y1*y1)
  663. r2=sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))
  664. if(r2<rup.AND.r2>=r)THEN
  665. Intensity=Intensity+s(x,y)
  666. zztop=zztop+1
  667. endif
  668. end do
  669. end do
  670. Intensity=Intensity/zztop
  671. ! if(r>0)THEN
  672. ! Circ=2*pi*r
  673. ! Intensity=(Intensity/Circ)
  674. ! endif
  675. WRITE(k,'(F12.4)')Intensity
  676. end do
  677. close(i)
  678. close(k)
  679. end do
  680. END DO
  681. END SUBROUTINE Gauss
  682. !SUBROUTINE PARAMETER(mint,rawin,tevout,home)
  683. !
  684. !INTEGER ::val,vala,valb
  685. !
  686. !235 Print *,"|---------------------------------------------------------------------|"
  687. !Print *,"| welcome, please select an option from the list below |"
  688. !Print *,"|---------------------------------------------------------------------|"
  689. !Print *,"| 0 -- Program Constants |"
  690. !print *,"| 1 -- File Directories |"
  691. !print *,"| 2 -- |"
  692. !print *,"| 2 -- |"
  693. !print *,"| 2 -- |"
  694. !print *,"| 5 -- |"
  695. !print *,"| 6 -- Exit to main menu |"
  696. !print *,"|_____________________________________________________________________|"
  697. !Read (*,*) val
  698. !if(val.eq.0)THEN
  699. !234 Print *,"|---------------------------------------------------------------------|"
  700. ! Print *,"| What Constant would you like to fiddle with? |"
  701. ! Print *,"|---------------------------------------------------------------------|"
  702. ! Print *,"| 0 -- 'MINT' the minimum threshold |"
  703. ! print *,"| 1 -- 'ARRAY SIZE' The size of the imeage input |"
  704. ! print *,"| 2 -- NA |"
  705. ! print *,"| 2 -- NA |"
  706. ! print *,"| 2 -- NA |"
  707. ! print *,"| 5 -- NA |"
  708. ! print *,"| 6 -- Exit to main menu |"
  709. ! print *,"|_____________________________________________________________________|"
  710. ! Read (*,*) vala
  711. ! IF(vala.eq.0)THEN
  712. ! Print*,"The current Minimum threshold is: ",mint
  713. ! print*,"Please enter the desired minimum threshold"
  714. ! READ(*,*)mint
  715. ! print*,"your value has been changed"
  716. ! goto 243
  717. ! elseif(valva.eq.1)THEN
  718. ! print*,"The current array size is: ",ARRAY_SIZE
  719. ! print*,"Please indicate the desired image size: "
  720. ! print*,'NOTE: The array begins at zero and therefore array size is the number of pixels minus one.'
  721. ! READ(*,*)ARRAY_SIZE
  722. ! elseif(vala.eq.6)
  723. ! endif
  724. !else if(val.eq.1) THEN
  725. !233 Print *,"|---------------------------------------------------------------------|"
  726. !Print *,"| Directories are tricky this menu is a long time in the making. |"
  727. !Print *,"|---------------------------------------------------------------------|"
  728. !Print *,"| 0 -- Raw Data Input Files |"
  729. !print *,"| 1 -- Time Evolution Output |"
  730. !print *,"| 2 -- Gaussian parameter |"
  731. !print *,"| 3 -- Home Directory (The highest common directory) |"
  732. !print *,"| 4 -- |"
  733. !print *,"| 5 -- |"
  734. !print *,"| 6 -- Exit to main menu |"
  735. !print *,"|_____________________________________________________________________|"
  736. !Read (*,*) valb
  737. !if(valb.eq.0)THEN
  738. !elseif(valb.eq.1)THEN
  739. ! PRINT*,"Not sure what to say here: ",rawin
  740. !elseif(valb.eq.2)THEN
  741. ! PRINT*,"funny you think this text is easy to write: ",tevout
  742. !elseif(valb.eq.3)THEN
  743. ! PRINT*,"Current home directory is: "home
  744. !elseif(valb.eq.6)THEN
  745. !endif
  746. !elseif(val.eq.6)THEN
  747. !endif
  748. !That should send us back to main by default
  749. !
  750. !
  751. !
  752. !
  753. !
  754. !
  755. !
  756. !END SUBROUTINE PARAMETER
  757. SUBROUTINE GAUSS_SUBTRACT(TABLER,TABLEI)
  758. !SAME OLD VARIABLES, NEW SONG!
  759. integer :: i,j,l,m,z,nums,spots,t,mjd,ljd,anyvariable,max,qq
  760. CHARACTER(len=3) :: cnum
  761. CHARACTER(len=2) :: cnum2
  762. CHARACTER(len=1) :: cnum1
  763. CHARACTER(len=3) :: posx
  764. CHARACTER(len=3) :: posy
  765. CHARACTER(len=4) :: n3
  766. CHARACTER(len=33) :: n1,daftpunk
  767. CHARACTER(len=1) :: n6
  768. CHARACTER(len=1) :: n7
  769. CHARACTER(len=1) :: u1
  770. CHARACTER(len=44) :: mp
  771. CHARACTER(len=40) :: nw
  772. REAL, DIMENSION(0:299,0:299) :: S,G,G2,F
  773. INTEGER :: XMAX,YMAX,XMIN,YMIN,XII,YII,x0,y0,dim1,dim2,k,x,y,xpeak,ypeak,responsum
  774. REAL :: Sum,I2,mint, fxy,x1,y1,xpeak1,ypeak1,sigma,mew
  775. LOGICAL :: hortari
  776. print*,"Please enter number of spots you wish to regress."
  777. Read(*,*)max
  778. Print*,"Would you like your values input from HEADER.dat? Note: 1 for yes and 2 for no will suffice."
  779. Read*,responsum
  780. if (responsum.eq.1)THEN
  781. OPEN(623,FILE='HEADER.dat') ! LETS be versitile here let's construct this file in the following way
  782. hortari=.TRUE.
  783. elseif(responsum.eq.0)THEN
  784. continue
  785. hortari=.FALSE.
  786. else
  787. Print*,"invalid response, You are now forced to manually input the values."
  788. hortari=.FALSE.
  789. endif
  790. n3='.dat'
  791. n1='/home/ramalab/research/dj11/dj11f' !5+8+9+5+6=
  792. daftpunk='/home/ramalab/research/regressed/'!5+8+9+11
  793. u1='_'
  794. ! 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.
  795. mjd=0
  796. ljd=0
  797. do anyvariable=0,max
  798. qq=anyvariable
  799. !print*,"Please input t, xpeak, ypeak, amp, sigma, and mew"
  800. IF(hortari.eqv..TRUE.)THEN
  801. 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
  802. Print *,t,xpeak,ypeak,amp,sigma,mew
  803. else
  804. PRINT*," PLEASE input the values for time, xpeak, ypeak, Amplitude, Sigma and Mew manually now: "
  805. READ *,t,xpeak,ypeak,Amp,Sigma,Mew
  806. endif
  807. ! REAL, DIMENSION(0:dim1,0:dim2) :: V
  808. WRITE(posx,'(i3)')xmin
  809. WRITE(posy,'(i3)')ymin
  810. xmin=xpeak-5
  811. xmax=xpeak+5
  812. ymin=ypeak-5
  813. ymax=ypeak+5
  814. xpeak1=xpeak
  815. ypeak1=ypeak
  816. write(posx,'(i3)')xpeak
  817. write(posy,'(i3)')ypeak
  818. if(t<10)THEN
  819. WRITE(n6,'(i1)')mjd
  820. WRITE(n7,'(i1)')ljd
  821. write(cnum1,'(i1)')t
  822. nw=n1//n6//n7//cnum1//n3 !nw=33+1+1+1+4=40 n1=33
  823. mp=daftpunk//posx//u1//posy//u1//n6//n7//cnum1 !33+3+3+1+1+1+1+1=5+6+33=44
  824. else if(t>=10.AND.t<100)THEN
  825. WRITE(n6,'(i1)')mjd
  826. WRITE(cnum2,'(i2)')t
  827. nw=n1//n6//cnum2//n3
  828. mp=daftpunk//posx//u1//posy//u1//n6//cnum2
  829. else
  830. WRITE(cnum,'(i3)')t
  831. nw=n1//cnum//n3
  832. mp=daftpunk//posx//u1//posy//u1//cnum
  833. endif
  834. ! print *,nw
  835. OPEN(678,FILE=mp,STATUS='UNKNOWN',ACTION='WRITE')
  836. OPEN(qq,FILE=nw,STATUS='OLD',ACTION='READ')
  837. READ(qq,*)S
  838. !print*,xmin,ymin
  839. G=0
  840. do x=xmin,xmax
  841. do y=ymin,ymax
  842. !print*,"IN DO"
  843. x1=x
  844. y1=y
  845. r2=sqrt((x1*x1)+(y1*y1)) !not a vector
  846. r0=sqrt((xpeak1*xpeak1)+(ypeak1*ypeak1)) !still not a vector
  847. r=sqrt((x1-xpeak)*(x1-xpeak)+(y1-ypeak)*(y1-ypeak))
  848. f(x,y)=A*EXP(((r-mew)*(r-mew))/(sigma*sigma)) ! we can fix that later UBER PROBLEM! A- THESE AREN"T VECTORS
  849. G(x,y)=f(x,y)
  850. end do
  851. end do
  852. ! print *," OUT DO"
  853. ! S=S2 !why did I do this?
  854. !Print*,"DO THE DEW" ! we get this far befor ethe segfault?
  855. do x=0,299
  856. do y=0,299
  857. G2(x,y)=S(x,y)-G(x,y) ! we need original-G
  858. end do
  859. end do
  860. PRINT*,"Between Two loops with zack galifinakis"
  861. do x=0,299
  862. do y=0,299
  863. WRITE(678,'(F12.4)',ADVANCE='NO')G2(x,y)
  864. end do
  865. Write(678,*)
  866. end do
  867. !print*,"Dew did that do?"
  868. CLOSE(qq)
  869. CLOSE(678)
  870. end do
  871. !print*,"jewjitdew?"
  872. END SUBROUTINE GAUSS_SUBTRACT
  873. ! SUBROUTINE CONVERGE
  874. !END SUBROUTINE CONVERGE
  875. !This is written to be a possible subroutine added so that we can set up our alpha tables
  876. SUBROUTINE ARPHA(tabler,TABLEI)
  877. IMPLICIT NONE
  878. ! we need to find and index the peaks
  879. !these variables are for use with the table array and should not be used in any other array/loop
  880. REAL :: gnum,xpeak,ypeak,phipeak,rpeak
  881. REAL, DIMENSION(0:4,0:60) :: tableR
  882. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  883. INTEGER :: i,x,y,xmin,xmax,ymin,ymax,k
  884. !these are defined for everything else
  885. INTEGER p,qq,q,pq,qp,gs,g,gn,gg,ggp,gp,zz,SIZE,gngn,gr,gsgs,dim,LB,UB,BIN
  886. Real :: SUM, SUM2, SUM3, alphat, avgalpha,RISE,RUN,r0,x0,y0,INT,alphaerr,SIGMA,SIGMA1
  887. REAL,SAVE :: MINT=50.0
  888. REAL, DIMENSION(0:299,0:299) :: S
  889. REAL, DIMENSION(0:99) :: alpha
  890. CHARACTER(len=3) :: cnum
  891. CHARACTER(len=4) :: dat
  892. CHARACTER(len=33) :: n1
  893. CHARACTER(len=40) :: n2
  894. PRINT *,"Please enter the lower bound, upper bound, and number of bins: "
  895. READ(*,*)LB,UB,BIN
  896. k=0
  897. N1='/home/ramalab/research/dj11/dj11f'
  898. x0=142.0
  899. y0=169.0
  900. r0=sqrt(x0*x0+y0*y0)
  901. dat='.dat'
  902. ! MINT=80.0
  903. do i=0,60
  904. gnum=tablei(0,k)
  905. xpeak=tablei(1,k)
  906. ypeak=tablei(2,k)
  907. xmin=tablei(3,k)
  908. ymin=tablei(4,k)
  909. xmax=tablei(5,k)
  910. ymax=tablei(6,k)
  911. ! PRINT *,gnum,xpeak,ypeak
  912. gsgs=0
  913. gngn=BIN-1
  914. k=k+1
  915. do g=gsgs,gngn !this loop is strait forward to me before making changes I recommend sending me and email.
  916. gr=LB+5*g
  917. gg=gr+5
  918. ggp=UB-1*(gr-LB)
  919. gp=ggp-5
  920. SUM2=0
  921. SUM3=0
  922. ! PRINT *,"LOOP1 here I come"
  923. Loop1: Do p=gr,gg
  924. WRITE(CNUM,'(i3)')p
  925. n2=n1//cnum//dat !33+7 =40
  926. open(unit=23,file=n2,STATUS='OLD',ACTION='READ')
  927. READ(23,*)S
  928. dim=0
  929. SUM=0
  930. do y=ymin,ymax
  931. do x=xmin,xmax
  932. INT=S(X,Y)
  933. IF(INT>MINT)THEN
  934. SUM=SUM+S(x,y)
  935. DIM=DIM+1
  936. ENDIF
  937. end do
  938. end do
  939. SUM=SUM/DIM
  940. Sum2=SUM2+SUM
  941. close(23)
  942. ENDDO LOOP1
  943. ! PRINT *,"Loop 2 looking at you"
  944. Loop2: Do p=gp,ggp
  945. SUM=0
  946. WRITE(CNUM,'(i3)')p
  947. n2=n1//cnum//dat !33+7 =40
  948. open(unit=43,file=n2,STATUS='OLD',ACTION='READ')
  949. READ(43,*)S
  950. dim=0
  951. do y=ymin,ymax
  952. do x=xmin,xmax
  953. INT=S(X,Y)
  954. IF(INT>MINT)THEN
  955. ! print*,S(x,y)
  956. SUM=SUM+S(x,y)
  957. DIM=DIM+1
  958. ENDIF
  959. end do
  960. end do
  961. SUM=SUM/DIM
  962. Sum3=Sum3+sum
  963. close(43)
  964. ENDDO LOOP2
  965. ! PRINT *,"ZZ what you got for me"
  966. zz=g
  967. Rise=Sum3-Sum2
  968. Run=ggp-gg
  969. alpha(zz)=Rise/(5*Run)
  970. ! PRINT *,alpha(zz)
  971. enddo
  972. ! PRINT *,"end do now to"
  973. alphat=0
  974. do zz=0,BIN-1
  975. alphat=ALPHAT+alpha(zz)
  976. end do
  977. avgalpha=log(alphat)/log(BIN)
  978. ! print *,avgalpha
  979. TABLER(2,i)=avgalpha
  980. ! Error Calculation
  981. SIGMA1=0
  982. do zz=0,BIN-1
  983. SIGMA1=(log(alpha(zz))-log(avgalpha))*(log(alpha(zz))-log(avgalpha))+SIGMA1
  984. end do
  985. SIGMA=SIGMA1/BIN
  986. alphaerr=sqrt(SIGMA)
  987. TABLER(4,i)=alphaerr
  988. end do
  989. END SUBROUTINE ARPHA
  990. SUBROUTINE WRIT (TABLER,TABLEI)
  991. REAL, DIMENSION(0:4,0:60) :: TABLER
  992. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  993. INTEGER x,y
  994. OPEN(23,FILE='TABLE.dat',STATUS='NEW',ACTION='WRITE')
  995. do y=0,60
  996. do x=0,6
  997. WRITE(23,'(I4)', ADVANCE='NO')tableI(x,y)
  998. end do
  999. do x=0,4
  1000. WRITE(23,'(F10.4)',ADVANCE='NO')tableR(x,y)
  1001. end do
  1002. WRITE(23,*)
  1003. end do
  1004. CLOSE(23)
  1005. END SUBROUTINE WRIT
  1006. SUBROUTINE IMPORT(TABLER,TABLEI)
  1007. REAL, DIMENSION(0:4,0:60) :: TABLER
  1008. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  1009. INTEGER x,y,y5
  1010. OPEN(23,FILE='TABLE.dat',STATUS='OLD',ACTION='READ')
  1011. do y=0,60
  1012. 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)
  1013. READ(23,'(5F10.4)')tableR(0,y),tableR(1,y),tableR(2,y),tableR(3,y),tableR(4,y)
  1014. end do
  1015. CLOSE(23)
  1016. END SUBROUTINE IMPORT
  1017. SUBROUTINE SORT(TABLER,TABLEI)
  1018. REAL, DIMENSION(0:4,0:60) :: TABLER, TEMP
  1019. INTEGER, DIMENSION(0:6,0:60) :: TABLEI
  1020. INTEGER x,y
  1021. Real :: Phi0, dphi, Rpeak, Phipeak,r0,x1,y1,r1,r2,ref,radangle,minang,maxang,xmin,ymin,dphi1
  1022. Integer :: i, xpoint, ypoint, grp,y5, xmi, ymi
  1023. REAL,SAVE :: x0=142,y0=169
  1024. PRINT*,"PLEASE input the angle for spots to be sorted by, in radians."
  1025. READ(*,*)radangle
  1026. r0=sqrt(x0*x0+y0*y0)
  1027. do y=0,60
  1028. xpoint=TABLEI(1,y)
  1029. ypoint=TABLEI(2,y)
  1030. x1=xpoint
  1031. y1=ypoint
  1032. xmi=TABLEI(3,y)
  1033. ymi=TABLEI(6,y)
  1034. xmin=xmi
  1035. ymax=ymi
  1036. r1=sqrt(x1*x1+y1*y1)
  1037. r2=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
  1038. TABLER(1,y)=r2
  1039. phi=acos(abs(abs(x1-x0)/r2))
  1040. TABLER(0,y)=phi
  1041. dphi=abs(phi-acos(abs(x1-x0)/(sqrt((xmin-x0)*(xmin-x0)+(ymax-y0)*(ymax-y0)))))
  1042. TABLER(3,y)=dphi
  1043. end do
  1044. do y=0,40
  1045. phi0=TABLER(0,y)
  1046. grp=TABLEI(0,y)
  1047. dphi=TABLER(3,y)
  1048. do y5=0,40
  1049. phi1=TABLER(0,y5)
  1050. dphi1=TABLER(3,y5)
  1051. delta=abs(phi0-phi1)
  1052. maxang=radangle+abs(dphi/2)+abs(dphi1/2)
  1053. minang=radangle-abs(dphi/2)-abs(dphi1/2)
  1054. if(delta<maxang.AND.delta>minang)THEN
  1055. TABLEI(0,y5)=grp
  1056. else
  1057. endif
  1058. enddo
  1059. enddo
  1060. END SUBROUTINE SORT
  1061. SUBROUTINE Q_SPACE
  1062. IMPLICIT NONE
  1063. INTEGER,DIMENSION(0:6,0:60) :: TABLEI
  1064. REAL,DIMENSION(0:4,0:60) :: TABLER
  1065. INTEGER :: xmin,ymin,xmax,ymax,i,xpeak,ypeak,rup,rdown,r,j,spots,x,y,rmax,nums,k,zztop,circumd,circ,circum
  1066. REAL :: r0,Intensity,r1,r2,x0,y0,x1,y1,PI,dphi,phiup,phimax,phimin,phi,rmin
  1067. REAL, DIMENSION(0:299,0:299) :: s,q,p
  1068. ! Used to put variables in file names
  1069. CHARACTER(len=1) :: cnum ! 1-9
  1070. CHARACTER(len=2) :: cnum2 ! 10-99
  1071. CHARACTER(len=3) :: cnum3,PHIW !100-360
  1072. CHARACTER(len=4) :: posx !xvalue initially found at --\___might want to change
  1073. CHARACTER(len=4) :: posy !y value initially found at --/ to peak center
  1074. ! Used to mark directories
  1075. CHARACTER (len=33):: N1 !/HOME/MICHAEL/RESEARCH/DJ11/DJ11F+ OTHER
  1076. CHARACTER(len=34) :: nu
  1077. CHARACTER (len=38):: n2 !output path same as input+ other
  1078. CHARACTER (len=39):: n2a !output path same as input+ other
  1079. CHARACTER (len=40):: n2b !output path same as input+ other
  1080. CHARACTER(len=4) :: n3 !.dat
  1081. CHARACTER(len=1) :: n6 !used to add 0
  1082. CHARACTER(len=1) :: n7 ! used to add 0
  1083. CHARACTER(len=1) :: u1 ! used to add _
  1084. CHARACTER(len=55) :: n4 !Complete input path n1+ cnum and n3 +1 for _ added
  1085. CHARACTER (len=52):: n5 !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  1086. CHARACTER (len=51):: n5a !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  1087. CHARACTER (len=54):: n5b !Complete OUTPUT path n1+ cnum and n3 +1 for _ added
  1088. !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.
  1089. OPEN(unit=900,FILE='OUTPUT.data',STATUS='OLD',ACTION='READ')
  1090. OPEN(unit=901,FILE='1DGauss.data',STATUS='UNKNOWN',ACTION='WRITE') ! This will be the output file.
  1091. Print *,"Please indicate the number of spots in the file. NOTE: there is no direct input functionality at this time."
  1092. Read (*,*) nums
  1093. PRINT *,"Please indicate the maximum desired radius from the center of the peak"
  1094. READ (*,*)rmax
  1095. PRINT*,"Please indicate the number of phi bins to be used"
  1096. READ(*,*)circumd !from the latin for circle, circumdare
  1097. n3='.dat'
  1098. u1='_'
  1099. n1='/home/ramalab/research/dj11/dj11f'
  1100. nu='/home/ramalab/research/Q_spc/dj11Q'
  1101. PI=3.14159
  1102. x0=142
  1103. y0=169
  1104. Do spots=1,nums
  1105. Read (900,'(6I4)')xmin,ymin,xmax,ymax, xpeak, ypeak
  1106. write(posx,'(i4)')xpeak
  1107. write(posy,'(i4)')ypeak
  1108. 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
  1109. phimin=acos(sqrt((xmax-x0)*(xmax-x0)+(ymin-y0)*(ymin-y0))/x0)
  1110. phimax=acos(sqrt((xmin-x0)*(xmin-x0)+(ymax-y0)*(ymax-y0))/x0)
  1111. R0=sqrt((xpeak-x0)*(xpeak-x0)+(ypeak-y0)*(ypeak-y0))
  1112. ! LOOP 1 OF 3 !
  1113. circum=circumd-1
  1114. do i=1,9
  1115. !defining out coordinate system
  1116. do circ=0,circum
  1117. phiup=phimin+dphi*(circ+1)
  1118. k=i+400
  1119. WRITE(CNUM,'(i1)')i
  1120. WRITE(PHIW,'(F3.0)')phiup
  1121. n2=n1//cnum//n3 !33+5 =38
  1122. n5=nu//posx//posy//u1//cnum//u1//PHIW//n3 !34+7+1+4
  1123. OPEN(unit=k,file=n5,STATUS='new',ACTION='WRITE')
  1124. open(unit=i,file=n2,STATUS='OLD',ACTION='READ')
  1125. READ(i,*)s
  1126. do r=0,rmax
  1127. rmin=r
  1128. do x=xmin,xmax
  1129. do y=ymin,ymax
  1130. x1=x
  1131. y1=y
  1132. phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
  1133. r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
  1134. if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
  1135. Intensity=S(x,y)+Intensity
  1136. endif
  1137. end do
  1138. end do
  1139. Intensity=Intensity/2
  1140. ! if(r>0)THEN
  1141. ! Circ=2*PI*r
  1142. ! Intensity=Intensity/Circ
  1143. ! endif
  1144. WRITE(k,'(F12.4," ",I3)')Intensity,r
  1145. end do
  1146. CLOSE(i)
  1147. CLOSE(k)
  1148. enddo
  1149. enddo
  1150. !LOOP 2 OF 3
  1151. do i=10,99
  1152. do circ=0,circum
  1153. phiup=phimin+dphi*(circ+1)
  1154. k=i+400
  1155. WRITE(CNUM2,'(i2)')i
  1156. WRITE(PHIW,'(F3.0)')phiup
  1157. n2a=n1//cnum2//n3 !33+7 =40
  1158. n5a=nu//posx//posy//u1//cnum2//n3 !34+7 =44
  1159. OPEN(unit=k,file=n5a,STATUS='NEW',ACTION='WRITE')
  1160. open(unit=i,file=n2a,STATUS='OLD',ACTION='READ')
  1161. READ(i,*)s
  1162. do r=0,rmax
  1163. rmin=r
  1164. do x=xmin,xmax
  1165. do y=ymin,ymax
  1166. x1=x
  1167. y1=y
  1168. phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
  1169. r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
  1170. if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
  1171. Intensity=S(x,y)+Intensity
  1172. endif
  1173. end do
  1174. end do
  1175. Intensity=Intensity/2
  1176. ! if(r>0)THEN
  1177. ! Circ=2*PI*r
  1178. ! Intensity=Intensity/Circ
  1179. ! endif
  1180. WRITE(k,'(F12.4," ",I3)')Intensity,r
  1181. end do
  1182. CLOSE(i)
  1183. CLOSE(k)
  1184. enddo
  1185. enddo
  1186. ! LOOP 3 OF 3
  1187. do i=100,360
  1188. !defining out coordinate system
  1189. do circ=0,circum
  1190. phiup=phimin+dphi*(circ+1)
  1191. k=i+400
  1192. WRITE(CNUM3,'(i3)')i
  1193. WRITE(PHIW,'(F3.0)')phiup
  1194. n2b=n1//cnum3//n3 !33+7 =40
  1195. n5b=nu//posx//posy//u1//cnum3//u1//phiw//n3 !34+7 =44
  1196. open(unit=i,file=n2b,STATUS='OLD',ACTION='READ')
  1197. OPEN(unit=k,file=n5b,STATUS='NEW',ACTION='WRITE')
  1198. do r=0,rmax
  1199. rmin=r
  1200. do x=xmin,xmax
  1201. do y=ymin,ymax
  1202. x1=x
  1203. y1=y
  1204. phi=acos(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)))
  1205. r1=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
  1206. if(phi<phiup.AND.phi>phimin.and.r1<rmax.and.r1>=rmin)THEN
  1207. Intensity=S(x,y)+Intensity
  1208. endif
  1209. end do
  1210. end do
  1211. Intensity=Intensity/2
  1212. ! if(r>0)THEN
  1213. ! Circ=2*PI*r
  1214. ! Intensity=Intensity/Circ
  1215. ! endif
  1216. WRITE(k,'(F12.4," ",I3)')Intensity,r
  1217. end do
  1218. CLOSE(i)
  1219. CLOSE(k)
  1220. enddo
  1221. ENDDO
  1222. enddo
  1223. END SUBROUTINE Q_SPACE

comments powered by Disqus