-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathid31offsets.f90
2768 lines (2735 loc) · 101 KB
/
id31offsets.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module specfiles
integer(kind=4) :: iunit,iscan,ispecerr,ncolumns
integer(kind=4),parameter :: LINELENGTH=512, WORDLENGTH=100
integer(kind=4),parameter :: NWORDS=60, NLINES=50
integer(kind=4),dimension(NLINES) :: headerwords
real(kind=8) :: scanstart, scanend
real(kind=8),dimension(NWORDS,NLINES) :: wordvalues
real(kind=8),dimension(NWORDS) :: Q
character(len=6) :: lnfmt
character(len=LINELENGTH) :: filnam, line
character(len=LINELENGTH) :: specsfilename, specsdate
character(len=WORDLENGTH) :: scantype
character(len=WORDLENGTH),dimension(NWORDS,NLINES) :: words
character(len=WORDLENGTH),dimension(NWORDS) :: columnlabels
integer(kind=4) :: i0, i9, id, isp, im
character(len=WORDLENGTH) :: FIRSTDET,LASTDET,TWOTTH
character(len=WORDLENGTH) :: MONITORCOL="Monitor"
logical topofscan
data lnfmt /'(a512)'/
data iunit /15/
contains
subroutine getfile
logical :: od
! See if there is already a file open, if there is then close it
inquire(unit=iunit,opened=od)
if(od)close(iunit)
! filnam must already be filled in, or we'll try to open garbage
1 open(unit=iunit,file=filnam,form='formatted', &
& access='sequential', status='old',iostat=ispecerr)
if(ispecerr.ne.0) then
write(*,'(3a,i6)')'Problem with file ',filnam(1:len_trim(filnam)),&
& ' iostat=',ispecerr
write(*,'(a)')'Please supply a valid SPEC file name'
read(*,lnfmt)filnam
goto 1
endif
topofscan=.false.
! Supply defaults in case dumb specfile has no #L
ncolumns=14 ; columnlabels(1)='2_theta'; columnlabels(2)='Epoch'
columnlabels(3)='Seconds'; columnlabels(4)='MA0'
columnlabels(5)='MA1'; columnlabels(6)='MA2'
columnlabels(7)='MA3'; columnlabels(8)='MA4'
columnlabels(9)='MA5'; columnlabels(10)='MA6'
columnlabels(11)='MA7'; columnlabels(12)='MA8'
columnlabels(13)='Monitor'; columnlabels(14)='Fluo det'
i0=ichar('0');i9=ichar('9')
id=ichar('.');isp=ichar(' ');im=ichar('-')
return
end subroutine getfile
subroutine findscan(n)
integer(kind=4),intent(in)::n
character(len=7) :: rd ! enough space for yes, no and unknown
inquire(unit=iunit,read=rd) ! Can we read the file?
if(rd(1:3).ne.'YES') call getfile ! If not go open it
if(iscan.eq.n .and. topofscan) return
if(iscan.gt.n) rewind(n) !
2 call readheader ! in case top of file
if(iscan.ne.n)then ! work through file to find #
1 read(iunit,lnfmt,err=10,end=10)line
if(line(1:1).eq.'#') goto 2
go to 1
else
! Scan found - check scantype elsewhere - not a job for specfile module
return
endif
10 write(*,*)
write(*,'(a,i5,a)')'Scan ',n,' not found in file '// &
& filnam(1:len_trim(filnam))
ispecerr=-2
return
end subroutine findscan
subroutine getdata(a,n)
! A is an array, dimensioned by the number of data items to be expected
integer(kind=4),intent(in) :: n
real(kind=8),dimension(n),intent(out) :: a
integer(kind=4), parameter :: one = 1
if(topofscan)then; topofscan=.false. ; else
read(iunit,lnfmt,err=10,end=10)line ! goes via line
topofscan=.false.
endif
call rdnums(one,n,a) ! reads from line
if(line(1:1).eq.'#')then; ispecerr=-1 ;goto 100 ; endif
go to 100
10 ispecerr=-1 ! End of file
100 return ; end subroutine getdata
subroutine readheader
character(len=1) :: letter
character(len=128) :: motor
integer(kind=4) :: i,j
real(kind=8) :: a
integer(kind=4),parameter :: four=4
if(line(1:1).eq.'#') goto 2 ! if already on header, don't skip
1 read(iunit,lnfmt,end=100)line
2 if(line(1:1).eq.'#') then; letter=line(2:2); topofscan=.true.
select case((letter))
case('O')
read(line(3:3),'(i1)')i
call split(line(4:LINELENGTH),words(:,i+1),WORDLENGTH,NWORDS,j)
headerwords(i+1)=j
case('P')
read(line(3:3),'(i1)')i
call rdnums(four,NWORDS,wordvalues(:,i+1))
case('Q')
! Copy the Q from the header into our specfile module
read(line(3:LINELENGTH),*,end=1, err=1) Q
case('N')
read(line(3:len_trim(line)),*,end=1)ncolumns
case('L') ! signals end of header, bug out here !
call split(line(3:LINELENGTH),columnlabels,WORDLENGTH,NWORDS,i)
if(ncolumns.ne.i) &
& write(*,'(a,i5)')'error reading header for scan',iscan
case('S')
read(line(3:len_trim(line)),*,end=1)iscan,scantype,motor, &
& scanstart, scanend
case('F')
specsfilename=line(3:len_trim(line)) ! get original filename
case('D')
read(line(3:len_trim(line)),'(a256)',end=1)specsdate
case('C'); continue ! comments
case('G'); continue ! no idea - always zero
case('E'); continue ! epoch - we don't care what it was for now.
case default
end select
else ! Not a # line, so assume end of header
read(line,*,err=1,end=1)a ! catch blank lines in header
topofscan=.true. !!! Must be able to read a number !
return !!! escapes from routine here !!!!!
endif
goto 1
100 ispecerr=-1; return; end subroutine readheader
real(kind=8) function getheadervalue(string)
character(len=*),intent(in) :: string
integer(kind=4) :: i, j, k
k=len(string); getheadervalue=0.0
do i=1,NWORDS; do j=1,NLINES
if(string(1:k).eq.words(i,j)(1:k))then
getheadervalue=wordvalues(i,j); return
endif
enddo; enddo
return; end function getheadervalue
integer(kind=4) function whichcolumn(string)
! Interprets the #L line information
character(len=*),intent(in) :: string
integer(kind=4) :: i, j
j=len(string) ; whichcolumn=-1
do i=1,NWORDS
if(string(1:j).eq.columnlabels(i)(1:j))then
whichcolumn=i ; return
endif
enddo
return; end function whichcolumn
subroutine rdnums(ic,n,values)
! Placed here in a subroutine in case of formatting or error handling problems
integer(kind=4),intent(in) :: n, ic
real(kind=8),dimension(n),intent(inout) :: values
if(line(ic:ic).eq.'#')then; ispecerr=-1; return; endif
read(line(ic:len_trim(line)),*,err=10,end=20)values(1:ncolumns)
goto 100
10 write(*,*)'Error reading line, looking for ',ncolumns,' values'
write(*,*)line(ic:len_trim(line))
write(*,*)values
20 continue
100 return
end subroutine rdnums
subroutine split(instring,outstrings,lenout,n,i)
integer(kind=4),intent(in) :: lenout,n
character(len=*),intent(in) :: instring
character(len=lenout),dimension(n),intent(out) :: outstrings
integer(kind=4),intent(out) :: i
integer(kind=4) :: j,k,l
j=1; k=1
do i=1,len_trim(instring) ! hope len > len_trim or array overstep
if(instring(i:i+1).ne.' ')then
outstrings(j)(k:k)=instring(i:i)
if(k.ne.1 .and. k.lt.lenout) k=k+1
if(instring(i:i).ne.' '.and.k.eq.1)k=k+1
else
if(k.gt.1)then ; do l=k,lenout
outstrings(j)(l:l)=' ' ! blank pad end of string
enddo ; j=j+1; k=1; if(j.gt.n)exit ; endif
endif
enddo
! Blank pad last string if necessary
if(k.gt.1)then
do l=k,lenout;outstrings(j)(l:l)=' ';enddo
endif
i=j; return; end subroutine split
end module specfiles
module rebin
integer(kind=4), parameter :: NCHAN = 9 ! ID31 will have nine channels
integer(kind=4) :: NCHANNEL = 9 ! ID31 will have nine channels
real(kind=8), dimension(NCHAN) :: offset, mult, multerr
logical :: tempres
integer(kind=4)logexdet(NCHAN)
integer(kind=4) :: iexrc ! region count
integer(kind=4),allocatable, dimension(:) :: iexarray
real(kind=8),allocatable, dimension(:,:) :: exarray
data logexdet /9*0/ ! whole array initialised to false
data iexrc /0/ ! no excluded regions by default
real(kind=8), allocatable :: ascan(:,:)
real(kind=8) :: step, tthlow, tthhigh, aminstep
! if requested by the user
real(kind=8) :: user_step = 0.003d0
real(kind=8) :: user_tthlow = -30.0d0
real(kind=8) :: user_tthhigh = 160.0d0
integer(kind=4) :: npts
data tthlow, tthhigh, step, aminstep /-30.0d0, 160.0d0, 0.003d0, &
& 0.0002d0 /
real(kind=8) :: sumtotal(nchan), sumtotalmon, minmon=1.0d0
real(kind=8) :: winhigh,winlow,winhighread,winlowread
integer(kind=4) :: wincol
real(kind=8) :: minrenormsig=5.0
real(kind=8) :: randomstart = 0, randomval=0
real(kind=8) :: randomend = 0, randomvalend=0
character(len=90) :: wincnt
logical :: winlog=.false.
logical :: userandomstart = .false.
logical :: rstchan( nchan ) = .false.
logical :: userandomend = .false.
logical :: renchan( nchan ) = .false.
logical :: wavelength_set = .false.
real(kind=8) :: wavelength=0.0d0
! T = Two theta
! Q = 4 * pi * sin( two_theta * pi / 360.0 ) / wavelength
! R = [q squared] = Q * Q
character(len=1) :: units = 'T'
real(kind=8) :: four_pi, pi_over_360
contains
integer(kind=4) function ibin(tth)
real(kind=8), intent(in) :: tth
real(kind=8)::x
if(tth.ge.tthlow .and. tth.le.tthhigh) then
x=(tth-tthlow)/step+0.5d0
ibin=int(x)
else ; ibin=-1 ; endif ! Out of range for this ascan array
return ; end function ibin
real(kind=8) function tthhb(n)
integer(kind=4), intent(in) :: n
integer(kind=4)nplusone
nplusone=n+1
tthhb=tthlb(nplusone)
return; end function tthhb
real(kind=8) function tthlb(n)
integer(kind=4), intent(in) :: n
tthlb = (real(n,8)-0.5d0)*step + tthlow
return; end function tthlb
real(kind=8) function bincen(n)
integer(kind=4), intent(in) :: n
bincen = real(n,8)*step + tthlow
return; end function bincen
subroutine checkrebinpars
real(kind=8) :: x
if ((units.ne.'T').and. (.not. wavelength_set) ) return
if(step.lt.aminstep .and. (units.eq.'T'))then
write(*,'(a,G12.4)')'Step is a bit small, resetting to ',aminstep
step=aminstep
user_step = step
endif
x=tthlow/step
tthlow=real(int(x),8)*step
x=(tthhigh-tthlow)/step
npts=int(x)
tthhigh=tthhb(npts)
if (units .eq. 'T') then
write(*,'(3(G12.5,a),G12.5)')tthlow,' < tth < ',tthhigh, &
&' step=',step,' npts=',npts
endif
if (units.eq.'Q') then
write(*,'(3(G12.5,a),G12.5)')tthlow,' < 2pi/d < ',tthhigh, &
&' step=',step,' npts=',npts
endif
if (units.eq.'R') then
write(*,'(3(G12.5,a),G12.5)')tthlow,' < Q^2 < ',tthhigh, &
&' step=',step,' npts=',npts
endif
return
end subroutine checkrebinpars
subroutine initialiserebin
integer(kind=4) :: ierr, i
real(kind=8) :: junk
! real :: time1, time2
! constants to machine precision
four_pi = 8.0d0*asin(1.0d0)
pi_over_360 = asin(1.0d0)/180.0d0
! default values
! Germanium numbers
offset(1) = 7.88643510d0
offset(2) = 5.91013889d0
offset(3) = 3.89405184d0
offset(4) = 1.97036951d0
offset(5) = 0.00000000d0
offset(6) = -2.12832415d0
offset(7) = -4.03585040d0
offset(8) = -6.00076222d0
offset(9) = -8.03007953d0
! Silicons latest numbers
offset(1) = 8.0634515d0
offset(2) = 5.8865159d0
offset(3) = 3.9594961d0
offset(4) = 2.0986688d0
offset(5) = 0.0000000d0
offset(6) = -1.9488783d0
offset(7) = -3.9966086d0
offset(8) = -6.0474594d0
offset(9) = -8.0536348d0
! Some new numbers for silicon - there seems to be some drift?
offset(1) = 8.05624406d0
offset(2) = 5.88332391d0
offset(3) = 3.95657893d0
offset(4) = 2.09530059d0
offset(5) = 0.00000000d0
offset(6) = -1.94883199d0
offset(7) = -3.99982480d0
offset(8) = -6.04725698d0
offset(9) = -8.05585677d0
! some more new numbers for silicon - there is drift for sure
offset(1) = 8.02703050d0
offset(2) = 5.88348041d0
offset(3) = 3.95722668d0
offset(4) = 2.09585757d0
offset(5) = 0.00000000d0
offset(6) = -1.94681946d0
offset(7) = -3.99878112d0
offset(8) = -6.04566287d0
offset(9) = -8.05515342d0
open(unit=16,file='temp.res',form='FORMATTED', &
& access='SEQUENTIAL', status='OLD', iostat=ierr)
if(ierr.eq.0)then
do i=1,nchannel
! Should clarify policy on whether to read/use these?
read(16,*,err=10,end=12)offset(i),mult(i),multerr(i),junk
enddo
! Report temp.res found and read in
write(*,'(a)')'temp.res file found and read in'
tempres=.true.
goto 11
endif
10 write(*,'(a)')'Not able to read temp.res file,'
write(*,'(a)')'You need this file for the detector calibration'
write(*,'(a)')'Please copy this file to the current directory'
write(*,'(a)')' eg: " cp ~/temp.res ." '
stop
12 write(*,'(a)')'Reached end of temp.res file early?'
11 close(16)
! two theta low, two theta high, step and that npts is correct
if ((units.eq.'T') .or. wavelength_set) call unit_lims
return ; end subroutine initialiserebin
subroutine rebinallocate
implicit none
integer ierr
if( allocated(ascan) ) deallocate(ascan) ! can reset
allocate(ascan(nchannel*2,npts),stat=ierr)
if(ierr.ne.0)stop 'Memory allocation error'
! call cpu_time(time1)
ascan=0.0d0 ! Clear any junk
! call cpu_time(time2)
! write(*,*)'Time taken to zero array =',time2-time1,'/s'
return ; end subroutine rebinallocate
subroutine bin(tth1,tth2,cts,ichan)
! tthh & tthl are the high & low two thetas from the SPEC file
! cts is the number of counts and ichan is the column to use in ascan
real(kind=8), intent(in) :: tth1, tth2, cts
integer(kind=4), intent(in) :: ichan
integer(kind=4) :: ibl, ibh, j
real(kind=8) :: frac, tthh, tthl
! tthl=min(tth1,tth2) ! Ensures ascending data
! tthh=max(tth1,tth2) ! strange bug on Fe3O4 exp??
if(tth1 .ge. tth2)then
tthl=tth2; tthh=tth1
else
tthl=tth1; tthh=tth2
endif
ibl = ibin(tthl) ! Bin of lower point
ibh = ibin(tthh) ! Bin of upper point
if ((ibl.le.0).or.(ibh.lt.ibl).or.(ibh.gt.npts))return ! error
if ( ibh .eq. ibl) then ! All in one bin
ascan(ichan,ibl) = ascan(ichan,ibl) + cts
return
elseif ( ibh .eq. ibl+1) then ! Spread over two bins
frac=(tthhb(ibl)-tthl)/(tthh-tthl) ! Fraction in bin 1
ascan(ichan,ibl)=ascan(ichan,ibl) + cts*frac
ascan(ichan,ibh)=ascan(ichan,ibh) + cts*(1.0d0-frac)
return
else
frac=(tthhb(ibl)-tthl)/(tthh-tthl) ! Fraction in bin 1
ascan(ichan,ibl)=ascan(ichan,ibl) + cts*frac
frac=(tthh-tthlb(ibh))/(tthh-tthl) ! Fraction in last bin
ascan(ichan,ibh)=ascan(ichan,ibh) + cts*frac
frac= step/(tthh-tthl) ! Fraction in middle bins
do j = ibl+1,ibh-1,1
ascan(ichan,j)=ascan(ichan,j) + cts*frac
enddo
return
endif
end subroutine bin
subroutine processline(tth1,tth2,cts,n,ctsmon)
real(kind=8), intent(in) :: tth1, tth2, ctsmon
integer(kind=4), intent(in) :: n
real(kind=8), intent(in), dimension(n) :: cts
integer(kind=4) :: i, ich
real(kind=8) :: tthh, tthl
do i = 1,n
if(logexdet(i).eq.1)cycle ! skip if excluded completely
if(userandomstart.and.rstchan(i).and.(tth2.lt.randomval))cycle
if(userandomend.and.renchan(i).and.(tth2.gt.randomvalend))cycle
tthl = tth1 - offset(i)
tthh = tth2 - offset(i)
if(logexdet(i).eq.2)then
! Is this an excluded region for this channel
if(led(i,tthl,tthh))cycle
endif
tthl = convert_unit_function( tthl )
tthh = convert_unit_function( tthh )
ich = 2*i-1
call bin(tthl,tthh,cts(i),ich)
ich = 2*i
call bin(tthl,tthh,ctsmon,ich)
sumtotal(i)=sumtotal(i)+cts(i)
enddo
sumtotalmon=sumtotalmon+ctsmon
return
end subroutine processline
logical function led(ic,xl,xh)
use specfiles
integer(kind=4),intent(in) :: ic
real(kind=8),intent(in) :: xl, xh
integer(kind=4)i
if(.not.allocated(iexarray).or. .not. allocated(exarray)) then
! oh dear, function should never have been called
write(*,'(a)')'Bug in program, bailing out, sorry!'
write(*,'(a)')'Please mail a bug report to [email protected]'
stop
endif
led=.false.
do i=1,iexrc
if((iexarray(i)+1).eq.ic)then ! found the right line
if(int(exarray(i,3)).lt.1)then
if(xl.gt.exarray(i,1) .and.xl.lt.exarray(i,2))led=.true.
if(xh.gt.exarray(i,1) .and.xh.lt.exarray(i,2))led=.true.
else
if(xl.gt.exarray(i,1) .and.xl.lt.exarray(i,2) .and. &
& iscan.eq.int(exarray(i,3)))led=.true.
if(xh.gt.exarray(i,1) .and.xh.lt.exarray(i,2) .and. &
& iscan.eq.int(exarray(i,3)))led=.true.
endif
endif
enddo
return
end function led
logical function window(ar,n)
use specfiles !
integer(kind=4),intent(in) :: n
real(kind=8),dimension(n),intent(in) :: ar
if(winlog)then
wincol=whichcolumn(wincnt(1:len_trim(wincnt)))
if(wincol.lt.0)then
write(*,*)'PROBLEMS WITH YOUR WINDOW - cannot find column', &
& wincnt(1:len_trim(wincnt))
write(*,*)'No further windowing will be attempted'
winlog=.false.
window=.true.
else
if( (ar(wincol).gt.winlow) .and. (ar(wincol).lt.winhigh)) then
window=.true.
else
window=.false.
endif
endif
else
window=.true.
endif
end function window
real(kind=8) function convert_unit_function( x )
real(kind=8), intent(in) :: x
real(kind=8) :: qt
select case((units))
case('T') ! Two theta, do nothing
convert_unit_function = x
return
case('Q') ! Q - constants from initialise_rebin
convert_unit_function = four_pi*sin(x*pi_over_360)/wavelength
return
case('R') ! Q^2
qt = four_pi * sin( x * pi_over_360 )/wavelength
! use a signed quantity
! write(*,*)qt,qt*qt,SIGN( qt * qt , qt ),x
convert_unit_function = sign( qt * qt , qt )
return
case default
convert_unit_function = x
return
end select
return
end function convert_unit_function
!
!
! To be called once per scan - picks up wavelength from command line
subroutine convert_unit_setupQ( Q )
! from rebin module :::: logical :: wavelength_set = .false.
real(kind=8), dimension(60), intent(in) :: Q
if ( (units .eq. 'T') .or. wavelength_set) then
return
endif
if( Q(4) .gt. 1.0E-15 ) then
wavelength = Q(4)
! Assume always the same for all scans:
write(*,*)
write(*,*)'Got wavelength',wavelength,'from spec file Q[3]'
wavelength_set = .true.
endif
! Check for errors
! let them have Angstroms as meters (1e-11)
if ((wavelength .lt. 1.0E-15) .and. (units .ne. 'T') ) then
write(*,*)'Your wavelength is a bit small ',wavelength
write(*,*)'Giving up, try putting wvln=1.234 on command line'
stop
endif
! write(*,*)'calling setupQ from convert_unit_setupQ'
call unit_lims()
return
end subroutine convert_unit_setupQ
subroutine unit_lims()
real(kind=8) :: stmp
! apply the limits according to if the user changed them
tthlow = convert_unit_function( user_tthlow )
tthhigh = convert_unit_function( user_tthhigh )
! set the step size to be right at 30 degrees twotheta
step = convert_unit_function( 30.0d0+user_step )
step = step - convert_unit_function( 30.0d0 )
! round it to be some printable represention
write(*,*)
! round to something nicely printable
if(units.ne.'T')then
stmp = int(log10(step)-2.0d0)*1.0d0
stmp = exp( log(10.0d0) * stmp )
step = int(step/stmp)*stmp
endif
! write(*,*)'from',tthlow,'to',tthhigh,'step',step
call checkrebinpars
call rebinallocate
return
end subroutine unit_lims
real(kind=4) function rlcg( )
! This returns a float in the range 0 to 1 which is a
! deterministic sequence. Will be the same over many runs.
integer(kind=4) :: a = 69069, c = 1327217885
integer(kind=4), save :: x = 1
real(kind=4), save :: m = 2147483647
x = x*a + c
rlcg = (real(x)/m+1.0)/2.0
! write(*,'(I12,1X,F12.10)')x, rlcg
end function rlcg
end module rebin
module summation
integer(kind=4):: isc=0, i6s=0 ! which scale factor, 6sig pts
integer(kind=4)::nzap=0, nsuperzap=0 ! for zapping
real(kind=8),allocatable :: sumdata(:,:),hist(:,:)
real(kind=8) :: alp=0.5d0 ! Bayesian zero counts fudge
real(kind=8) :: scalinp=1.0d5 ! Scale factor for .inp files
real(kind=8) :: zap=6.0d0 ! level for outlier elimination (median filter?)
real(kind=8) :: superzaplevel=1.0d0 ! level for zinger elimination
logical :: renorm=.false. ! to update temp.res file efficiencies
logical :: zapping=.false.
logical :: superzap=.false.
logical :: medianofchannels=.false.
contains
subroutine calibsum
use specfiles
use rebin
integer(kind=4) :: m,n,ierr
real(kind=8),dimension(NCHAN) :: tmult,tmulterr
call effic(n,m,tmult,tmulterr)
call reporteffic(n,m,tmult,tmulterr)
ierr=0
if(.not.allocated(sumdata))allocate(sumdata(3,npts),stat=ierr)
! 3 ! cts, mon, e(mon)
if(ierr.ne.0)then
write(*,'(a)') 'Memory allocation error in calibsum'
stop
endif
if(.not.allocated(hist))allocate(hist(2,npts),stat=ierr)
! final signal and esd
if(ierr.ne.0)then
write(*,'(a)') 'Memory allocation error in calibsum'
stop
endif
call sumthem ! combine detectors in sumdata array
if(medianofchannels)call medianchannels
call normerr ! determine error bars and fill in hist
call checkdets ! check ascan matches hist
1 if(zapping)then
write(*,*)
write(*,'(A,F8.5,A)')"After zapping a sigma level ",zap," ... "
call sumthem ! combine detectors in sumdata array
call normerr ! determine error bars and fill in hist
call checkdets ! check ascan matches hist
endif
if(nzap.gt.0)goto 1
2 if(superzap .and. nsuperzap.gt.0)then
write(*,*)
write(*,'(A,F8.5,A)')"After superzapping at level ", superzaplevel,"..."
call sumthem ! combine detectors in sumdata array
call normerr ! determine error bars and fill in hist
call checkdets ! check ascan matches hist
endif
if(nsuperzap.gt.0)goto 2
return
end subroutine calibsum
subroutine effic(n,m,tmult,tmulterr)
use rebin ! mult, multerr, nchan, ascan, tempres, npts
integer(kind=4),intent(inout)::n,m
real(kind=8),intent(inout),dimension(NCHAN)::tmult,tmulterr
integer(kind=4) :: i,j,k,nex
real(kind=8), dimension(4,NCHAN) :: signal
real(kind=8) :: sumsig,sumsige
signal=0.0d0; n=0
write(*,'(a)')'Determining detector efficiencies'
nex=0; do j=1,nchannel; if(logexdet(j).eq.1)nex=nex+1; enddo
do i=1,npts
k=0
do j=1,nchannel ! must have more than 1 monitor in bin to use
if(ascan(2*j,i).gt.1.0d0)then
if(ascan(2*j-1,i).gt.minrenormsig)k=k+1
endif
enddo
if(k.eq.(nchannel-nex))then
n=n+1
do j=1,nchannel ! signal is a big bin for all overlapping points
if(logexdet(j).eq.1)cycle
! If a channel is excluded for the whole tth range in a file then this will fail?
signal(1,j)=signal(1,j)+ascan(2*j-1,i) ! counts
signal(2,j)=signal(2,j)+ascan(2*j,i) ! monitors
enddo
endif
enddo
if(n.gt.0)then
do j=1,nchannel ! Normalise Signal and get the error bar on it
if(logexdet(j).eq.1)cycle
signal(3,j)=signal(1,j)/signal(2,j)
signal(4,j)=signal(3,j) * &
& sqrt(1.0d0/signal(1,j)+1.0d0/signal(2,j))
enddo
sumsig=sum(signal(3,:)) ! sum of all signals in overlap region
sumsige=0.0d0
do i=1,nchannel ! Error in sum of all signals
sumsige=sumsige+signal(4,i)*signal(4,i)
enddo
sumsige=sqrt(sumsige)
do j=1,nchannel ! Finally get the channel efficiencies
if(logexdet(j).eq.1)then
tmult(j)=1.0d0; tmulterr(j)=0.0d0
else
tmult(j)=(real((nchannel-nex),8)*signal(3,j)/sumsig)
tmulterr(j)=tmult(j)*sqrt((signal(4,j)/signal(3,j))**2 + &
& (sumsige/sumsig)**2)
endif
enddo
if(.not.tempres)then ! Copy these to rebin module if no temp.res file
mult=tmult
multerr=tmulterr
m=1
else ! tempres
m=2 ! Flag the need to print the temporary (unused vals)
endif ! tempres
else ! n.gt.0
m=1 ! no overlap so only one to print
if(.not.tempres)then
mult=1.0d0 ! Make sure defaults are always 1.0d0
multerr=0.0d0
endif
endif ! if n.gt.0
return
end subroutine effic
subroutine sumthem
use rebin ! for ascan and mult
integer(kind=4) :: i,j
sumdata=0.0d0 ! the big array where the sum of channels goes
do i=1,npts
do j=1,nchannel
! counts
sumdata(1,i)=sumdata(1,i)+ascan(2*j-1,i)
! normalised mon = m1*e1 + m2*e2 + ....
sumdata(2,i)=sumdata(2,i)+ascan(2*j,i)*mult(j)
! emon**2 =
sumdata(3,i)=sumdata(3,i)+ &
! c ( e**2 + c * <e>**2 )
& ascan(2*j,i) * ( mult(j)**2 + ascan(2*j,i)*multerr(j)**2)
enddo
! emon=sqrt(emon**2)
sumdata(3,i)=dsqrt(sumdata(3,i))
enddo
return
end subroutine sumthem
subroutine normerr
! alp, hist & sumdata available as member of summation
use rebin ! npts
real(kind=8):: msq, s
integer(kind=4) :: i, n, n3
n=0; s=0.0d0; n3=0
do i=1,npts
if(sumdata(2,i).gt.minmon)then ! needs at least 1. mon count
hist(1,i)=sumdata(1,i)/sumdata(2,i) ! correct
msq=(sumdata(2,i)*sumdata(2,i)) ! always gt zero
hist(2,i)=sqrt( &
& (sumdata(1,i)+alp)/msq + (sumdata(1,i)*sumdata(3,i)/msq)**2)
n=n+1
s=s+hist(1,i)**2/hist(2,i)**2
if(hist(1,i)/hist(2,i).lt.3.0d0)n3=n3+1
else
hist(1,i)= 0.0d0 ! unobserved regions are filled with nonsense
hist(2,i)=-1.0d0
endif
enddo
write(*,1000)100.0d0*sqrt(real(n,8)/s),n3,n
1000 format('R_exp = ',f7.3,' with ',i7, &
& ' pts having I/<I> less than 3, from ',i7,' pts obs')
end subroutine normerr
integer(kind=4) function ctchan(i)
use rebin
integer(kind=4),intent(in) :: i
integer(kind=4) :: j
ctchan=0
do j=1,nchannel
if(ascan(2*j,i).gt.minmon)ctchan=ctchan+1
enddo; return; end function ctchan
subroutine checkdets
use rebin ! for ascan and mult arrays
integer(kind=4) :: i, j, ipt3s, ipt6s, iapts, myctchan
real(kind=8) :: c,y,ey2,wd2 ! chi2,sig,error^2 and wtd diff^2
real(kind=8),dimension(nchan)::sz ! superzap array
ipt3s=0 ; ipt6s=0; c=0.0d0; iapts=0; nzap=0; nsuperzap=0;
do i = 1, npts
myctchan=ctchan(i)
if(myctchan.ge.2)then
if(superzap)sz=-1.0
do j = 1, nchannel
if(ascan(2*j,i).gt.1.0d0)then ! need one monitor count to bother
y=ascan(2*j-1,i)/ascan(2*j,i)/mult(j) ! OK
! assumes zero error in the efficiency
ey2= ((ascan(2*j-1,i)+alp)/ascan(2*j,i)**2 &
& + ascan(2*j-1,i)**2/ascan(2*j,i)**3)/mult(j)**2
wd2=(y-hist(1,i))**2/(ey2+hist(2,i)**2)
if(superzap)sz(j)=y !
c=c+wd2
iapts=iapts+1
if(wd2.gt.9.0d0) ipt3s=ipt3s+1 ! 3s^2=9.0
if(wd2.gt.36.0d0)ipt6s=ipt6s+1 ! 6s^2=36.0
if(wd2.gt.zap*zap .and. zapping)then
ascan(2*j,i)=0.0 ! set mon and det to zero for zapped points
ascan(2*j-1,i)=0.0
nzap=nzap+1
endif
endif
enddo
! now check with superzap
if(superzap .and. myctchan.ge.3)call superzapem(sz,ascan(:,i),nchannel)
endif
enddo
if(iapts.gt.0) then
c=c/real(iapts,8)
write(*,1000)c,iapts
1000 format('Reduced chi**2 for channel merge =',F9.4,' from ', &
& i10,' pairs of pts')
write(*,1001)100.0d0*real(ipt3s)/real(iapts), &
&100.0d0*real(ipt6s)/real(iapts)
1001 format(f5.2,'% differ by >3 sigma, ',f7.4, &
& '% by >6 sigma (ideally 0.04% and 0.0000%)' )
i6s=ipt6s ! module variable
if(zapping)then
write(*,'(A,I8,A)')'Zapping ',nzap,' points'
endif
if(superzap)write(*,'(A,I8,A)')'Superzapping ',nsuperzap,' points'
else
write(*,*)'No channel overlap, whatsoever, was found'
i6s=0
endif
return
end subroutine checkdets
subroutine reporteffic(n,m,tmult,tmulterr)
use rebin ! mult, multerr, NCHAN
integer(kind=4),intent(in) :: n,m
real(kind=8),intent(in),dimension(nchan)::tmult,tmulterr
integer(kind=4) :: i
logical :: trex
if(n.gt.0)then
write(*,'(a,i9,a)')'Channel efficiences found from ',n, &
& ' points where all detectors overlap'
if(m.eq.1)then
write(*,'(a)')'Det Offset Effic <Effic>'
do i=1,NCHANNEL
! i-1 instead of i to start at channel zero
write(*,'(i3,3(1X,F10.7))')i-1,offset(i),mult(i),multerr(i)
enddo
endif ! m.eq.1
if(m.eq.2)then
write(*,'(a)')'Efficiencies from temp.res file,'// &
& ' the values found now are compared'
if(.not.renorm)then ; write(*,'(a)') &
& 'Det Offset Effic <Effic> current unused values'
else ; write(*,'(a)') &
& 'Det Offset Old Eff Old <E> New Eff New <E>'
endif
do i=1,NCHANNEL
! i-1 instead of i to start at channel zero
write(*,'(i3,5(1X,F10.7))')i-1,offset(i),mult(i),multerr(i), &
& tmult(i),tmulterr(i)
enddo
if(renorm)then ; mult=tmult ; multerr=tmulterr ; endif
endif ! m.eq.2
inquire(file='temp.res',exist=trex)
if(.not.trex .or. renorm)call tempreswrite
if(renorm) renorm=.false. ! for sumall - can only renorm once
else ! n.gt.0
write(*,'(a,$)')'No detector overlap found, efficiencies'
if(tempres)then
write(*,'(a)')' from temp.res file'
else
write(*,'(a)')' are probably wrong'
call tempreswrite
endif
write(*,'(a)')'Det Offset Effic <Effic>'
do i=1,NCHANNEL
! i-1 instead of i to start at channel zero
write(*,'(i3,3(1X,F10.7))')i-1,offset(i),mult(i),multerr(i)
enddo
endif
return
end subroutine reporteffic
subroutine tempreswrite
use rebin ! offset, nchan, mult
integer(kind=4) :: ierr, i
open(unit=16,file='temp.res',status='UNKNOWN',access='SEQUENTIAL',&
& form='FORMATTED',iostat=ierr)
if(ierr.ne.0)then
write(*,'(a)')'Couldn''t open temp.res to write!!!'
return !! bugs out
endif
write(*,'(a)')'Created temp.res file'
do i=1,nchannel
write(16,1000)offset(i),mult(i),multerr(i),0.0d0
enddo
1000 format(4(f11.8,','))
close(16)
tempres=.true. ! Should exist and be readable now
return; end subroutine tempreswrite
subroutine scaltot
! uses hist and sumdata arrays in summation module
real(kind=8) :: factor, sum1, sum2
sum1=sum(sumdata(1,:)) ! sum up real counts
sum2=sum(hist(1,:)) ! sum up normed counts
factor=sum1/sum2 ! determine multiplier
hist=hist*factor ! rescale hist array
return
end subroutine scaltot
subroutine scalpk
! uses hist and sumdata arrays in summation module
real(kind=8) :: factor
integer(kind=4), dimension(3) :: i
i=maxloc(sumdata,dim=2) ! get index of highest number of counts
factor=sumdata(1,i(1))/hist(1,i(1)) ! get scaling factor
hist=hist*factor
return
end subroutine scalpk
subroutine superzapem(sz,as,n)
integer(kind=4), intent(in) :: n
real(kind=8),dimension(n),intent(in) :: sz ! signal
real(kind=8),dimension(2*n),intent(inout) :: as ! ascan array
integer(kind=4) :: i, ib, nc
real(kind=8) :: s, ss, mean1, mean2, esd1, esd2, pc
! make mean and esd from sz
nc=0;s=0.0;ss=0.0
ib=1; pc=0.0
do i = 1,n
if(sz(i).ge.0.0)then
s=s+sz(i)
ss=ss+sz(i)*sz(i)
nc=nc+1
if(sz(i).gt.sz(ib))ib=i
endif
enddo
mean1=s/nc
esd1=(ss-mean1*mean1)/nc
! now without ib
s=0.0;ss=0.0;nc=0
do i = 1,n
if(sz(i).ge.0.0 .and. i.ne.ib)then
s=s+sz(i)
ss=ss+sz(i)*sz(i)
nc=nc+1
endif
enddo
mean2=s/nc
esd2=(ss-mean2*mean2)/nc
if(esd1.gt.0.0) pc=(esd1-esd2)/esd1 ! should use mean info too.
if(pc.gt.superzaplevel)then ! zap the ib channel
as(2*ib)=0.0
as(2*ib-1)=0.0
nsuperzap=nsuperzap+1
endif
return
end subroutine superzapem
! netlib code
! *DECK DSORT
SUBROUTINE DISORT (DX, DY, N, KFLAG)
!C***BEGIN PROLOGUE DSORT
!C***PURPOSE Sort an array and optionally make the same interchanges in
!C an auxiliary array. The array may be sorted in increasing
!C or decreasing order. A slightly modified QUICKSORT
!C algorithm is used.
!C***LIBRARY SLATEC
!C***CATEGORY N6A2B
!C***TYPE DOUBLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
!C***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING
!C***AUTHOR Jones, R. E., (SNLA)
!C Wisniewski, J. A., (SNLA)
!C***DESCRIPTION
!C
!C DSORT sorts array DX and optionally makes the same interchanges in
!C array DY. The array DX may be sorted in increasing order or
!C decreasing order. A slightly modified quicksort algorithm is used.
!C
!C Description of Parameters
!C DX - array of values to be sorted (usually abscissas)
!C DY - array to be (optionally) carried along
!C N - number of values in array DX to be sorted
!C KFLAG - control parameter
!C = 2 means sort DX in increasing order and carry DY along.
!C = 1 means sort DX in increasing order (ignoring DY)
!C = -1 means sort DX in decreasing order (ignoring DY)
!C = -2 means sort DX in decreasing order and carry DY along.
!C
!C***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm
!C for sorting with minimal storage, Communications of
!C the ACM, 12, 3 (1969), pp. 185-187.
!C***ROUTINES CALLED XERMSG
!C***REVISION HISTORY (YYMMDD)
!C 761101 DATE WRITTEN
!C 761118 Modified to use the Singleton quicksort algorithm. (JAW)
!C 890531 Changed all specific intrinsics to generic. (WRB)
!C 890831 Modified array declarations. (WRB)
!C 891009 Removed unreferenced statement labels. (WRB)
!C 891024 Changed category. (WRB)
!C 891024 REVISION DATE from Version 3.2
!C 891214 Prologue converted to Version 4.0 format. (BAB)
!C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
!C 901012 Declared all variables; changed X,Y to DX,DY; changed
!C code to parallel SSORT. (M. McClain)
!C 920501 Reformatted the REFERENCES section. (DWL, WRB)
!C 920519 Clarified error messages. (DWL)
!C 920801 Declarations section rebuilt and code restructured to use
!C IF-THEN-ELSE-ENDIF. (RWC, WRB)
!C***END PROLOGUE DSORT
!C .. Scalar Arguments ..
INTEGER(kind=4) KFLAG
INTEGER(kind=4) N
!C .. Array Arguments ..
DOUBLE PRECISION DX(*)
INTEGER(kind=4) DY(*)
!C .. Local Scalars ..
DOUBLE PRECISION R, T, TT, TTY, TY
INTEGER I, IJ, J, K, KK, L, M, NN
!C .. Local Arrays ..
INTEGER IL(21), IU(21)
!C .. External Subroutines ..
!C EXTERNAL XERMSG
!C .. Intrinsic Functions ..
INTRINSIC ABS, INT
!C***FIRST EXECUTABLE STATEMENT DSORT
NN = N
IF (NN .LT. 1) THEN
! CALL XERMSG ('SLATEC', 'DSORT',
! + 'The number of values to be sorted is not positive.', 1, 1)
RETURN
ENDIF
!C
KK = ABS(KFLAG)
IF (KK.NE.1 .AND. KK.NE.2) THEN
! CALL XERMSG ('SLATEC', 'DSORT',
! + 'The sort control parameter, K, is not 2, 1, -1, or -2.', 2,
! + 1)
RETURN
ENDIF
!C