-
Notifications
You must be signed in to change notification settings - Fork 0
/
binit.tex
9378 lines (8943 loc) · 481 KB
/
binit.tex
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
\newcommand{\NWtarget}[2]{#2}
\newcommand{\NWlink}[2]{#2}
\newcommand{\NWtxtMacroDefBy}{Fragment defined by}
\newcommand{\NWtxtMacroRefIn}{Fragment referenced in}
\newcommand{\NWtxtMacroNoRef}{Fragment never referenced}
\newcommand{\NWtxtDefBy}{Defined by}
\newcommand{\NWtxtRefIn}{Referenced in}
\newcommand{\NWtxtNoRef}{Not referenced}
\newcommand{\NWtxtFileDefBy}{File defined by}
\newcommand{\NWtxtIdentsUsed}{Uses:}
\newcommand{\NWtxtIdentsNotUsed}{Never used}
\newcommand{\NWtxtIdentsDefed}{Defines:}
\newcommand{\NWsep}{${\diamond}$}
\newcommand{\NWnotglobal}{(not defined globally)}
\newcommand{\NWuseHyperlinks}{}
% Nuweb formatted latex file for ID31 data reduction routines
% Most of this is bog standard latex
% Anything to do with @ characters is specific to nuweb
%
%
% The word FIXME anywhere in this document indicates
% an area where more attention is still needed.
\documentclass[10pt,a4paper,notitlepage]{article}
\usepackage{graphics} % For the pictures
\usepackage{anysize} % Try to circumvent Latex default margins
\newcommand{\var}[1]{\textbf{\textsf{#1}}} % highlight variables in text
\newcommand{\code}[1]{\textbf{\textsf{#1}}} % highlight code in text
\newcommand{\param}[1]{\textbf{\textsf{#1}}} % ... parameters ...
\newcommand{\mod}[1]{\textbf{\textsf{#1}}} % ... module names ...
\begin{document}
\marginsize{2cm}{2cm}{2cm}{2cm} % Needs anysize
\pagestyle{headings} % These are ugly - fix them somehow?
\pagenumbering{arabic} % Page numbers
\title{Processing powder data}
\author{JPW}
\date{Spring 2002}
\maketitle
\begin{abstract}
The real aim of this project was for the current author
to have a go at literate programming. As a side effect it is hoped that a
documented set of data reduction routines for ID31 will be produced. These
routines should not only be correct, but also include documentation to
inform the user what they do and how they do it, in excruciating detail.
They should fully replace the programs used for BM16, and hopefully be
understandable and extensible enough to accomodate any future developments.
\end{abstract}
\tableofcontents
% \newpage % Only need this if contents don't run over first page
\section{Introduction}
Powder diffraction data were traditionally recorded by moving a detector
to a particular angle, recording the number of counts produced in a given
time, and then moving on to measure at the next angle of interest.
The resulting data are then a set of angular positions
and counts, and perhaps esds (or statistical weights).
Conventionally these data are then analysed
by programs which insist on constant angular step sizes.
Collecting data in this format by step scanning has a disadvantage,
the detector arm has to keep stopping and starting, which is fairly
time consuming.
Also, using more than one detector makes it very difficult to ensure
data from different channels will differ by an integer number of steps.
The method used at ESRF is to scan the diffractometer at constant speed
and record the counts arriving and angular positions as a
function of time and then to take that (large) dataset and process it
into a conventional format for the user.
A series of scans can also be combined together, allowing greater
flexibility in experimental design.
This document is intended to illuminte the mysteries of the data
processing carried out.
\section{Specifications}
The program(s) will need to do the following things:
\begin{itemize}
\item Put the data on a constant step sized two theta scale
\item Combine data from multiple scans (and SPEC files even)
\item Check scans for consistency ($R_{merge}$) when combining them together
\item Determine detector offsets
\item Determine detector efficiency corrections
\item Apply any other corrections, for example flat plate stuff
\item Produce two dimensional datasets versus external variables
(eg temperatures)
\item Identify and deal with problems (gliches, zingers, shclose etc etc)
\item Demonstrate some sort of correctness
\item Carry out the tasks efficiently (be fast)
\item Doubtless other things currently undreamed of... (be vaguely extensible)
\end{itemize}
The data arrives in files, called SPEC files, which contain motor positions,
counts, times and perhaps temperatures.
Getting from this information to the final constant step sized single
histogram involves some extra information.
The angular offset and efficiency of each detector need to be known or
determined.
The counts need to be rebinned onto a commensurate two theta scale
and corrected for detector efficiency and summed together.
Finally some error bars for the data points need to be determined and
an overall scaling can be applied to the histogram.
\section{Program structure}
Having one large program to do all of these jobs is likely to result
in somethingwhich is as fast as possible, but might lack flexibility
in the future.
The intention is to produce a series of smaller programs and subroutines
which will work together and can also be combined into monolithic beasts.
(FIXME) This document is initially going to be a bit of a mess, but the hope
is that by being forced to improve and think
about this as a readable account of how the binning is done,
the programs will benefit.
We can delete that sentence once we're finished!
Fortran 90 is being used throughout, with no apologies from the author.
Knowledge of that language would be handy for anyone modifying these routines
in the future. Ideally there will be enough text and comments in here to
make things fairly self explanatory.
The code should be standard conforming (at least as far as the compilers
can enforce the standard).
Porting to fortran 77 probably just means doing
something about the allocatable arrays and removing some syntax.
Adding some fortran 95 bits would be desirable if we can get hold of a
compiler that generates parallel code.
The code is in something of a condensed format in places, with multiple
statements on a single line.
Given the amount of explanation in the accompanying text it seemed
reasonable to save space. A global
find and replace of ; with carriage return plus six spaces would probably
fix this.
Real numbers are universally kind=8 (double precision) and integers are
single precision.
Results have been found to be numerically identical on both windows
and linux based systems, suggesting that any floating point problems
remain fairly subtle.
(FIXME) This section wants an overview of the programs from the users
perspective. Or maybe a user manual section for the programs in a userio
section?
We'll begin with a section about SPEC files and some routines to handle them.
Then move onto rebinning, describing what is actually
recorded by the instrument during an experiment and how to convert that
into something we are more used to. This will culminate in a program which can
take a single scan out of a SPEC file and write out binned files for each
of the detectors separately.
Once the data are on convenient 2$\theta$ scales they can be added
together (more easily).
Combining the data from different channels and different scans means just summing them together, with appropriate corrections for detector efficiencies
(and other nastiness).
The determination of the detector efficiency is logically carrired out here,
along with some checks for compatibility of the scans.
Some statistics, like $R_{merge}$ will be introduced along with the
calculation of the statistical weights for the final histograms.
Unfortunately the $R_{merge}$ calculation for a particular scan requires
you to already know how the total looks, so that needs access to both
individual scans and the summed total - can be done on two passes
of the SPEC file or by processing the scans separately and then together.
The only thing obviously missing is the calculation of the detector offsets.
Due to the design of the instrument these should never change - so this
gets a section to itself and a bit of a standalone method and program.
Peak fitting is so ubiquitous that some routines might be added to the
package to do this - mainly for strain scanning experiments, otherwise
it gets a bit too involved.
\section{Reading SPEC files}
Fortran code for doing this probably exists already somewhere at ESRF, but we
didn't find it.
Something is knocked together here which is specific to BM16/ID31 files.
It would probably make a lot more sense to have used the specfile
library which is designed for c programmers. Probably this could be wrapped
up somehow but it looked a bit daunting. For now, some
fortran specific stuff is used, but we might come back to doing it
properly later.
A cursory examination of the specfile library suggests it tries to index
all of the scans in the file before doing anything useful, which is
something we want to avoid.
When someone comes along with a 100MB file containing a thousand scans that
might take a long time, so we aim to traverse a specfile in a single pass,
collecting the information we want as we go along.
File IO is likely to be the time critical step in the whole operation although
some profiling should be done if things aren't fast to begin with.
\subsection{Some variables to hang on to}
Some shared information amongst the various specfile
routines is needed, so we'll make a list of what we want here.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap1}\raggedright\small
\NWtarget{nuweb3}{} $\langle\,${\it SPECVARS}\nobreak\ {\footnotesize {3}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ integer(kind=4) :: iunit,iscan,ispecerr,ncolumns@\\
\mbox{}\verb@ integer(kind=4),parameter :: LINELENGTH=512, WORDLENGTH=100@\\
\mbox{}\verb@ integer(kind=4),parameter :: NWORDS=60, NLINES=50@\\
\mbox{}\verb@ integer(kind=4),dimension(NLINES) :: headerwords@\\
\mbox{}\verb@ real(kind=8) :: scanstart, scanend@\\
\mbox{}\verb@ real(kind=8),dimension(NWORDS,NLINES) :: wordvalues@\\
\mbox{}\verb@ real(kind=8),dimension(NWORDS) :: Q@\\
\mbox{}\verb@ character(len=6) :: lnfmt@\\
\mbox{}\verb@ character(len=LINELENGTH) :: filnam, line@\\
\mbox{}\verb@ character(len=LINELENGTH) :: specsfilename, specsdate@\\
\mbox{}\verb@ character(len=WORDLENGTH) :: scantype@\\
\mbox{}\verb@ character(len=WORDLENGTH),dimension(NWORDS,NLINES) :: words@\\
\mbox{}\verb@ character(len=WORDLENGTH),dimension(NWORDS) :: columnlabels@\\
\mbox{}\verb@ integer(kind=4) :: i0, i9, id, isp, im@\\
\mbox{}\verb@ character(len=WORDLENGTH) :: FIRSTDET,LASTDET,TWOTTH@\\
\mbox{}\verb@ character(len=WORDLENGTH) :: MONITORCOL="Monitor"@\\
\mbox{}\verb@ logical topofscan@\\
\mbox{}\verb@ data lnfmt /'(a512)'/@\\
\mbox{}\verb@ data iunit /15/ @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
Those variables are intended to contain the following information: \var{iunit}
is the fortran unit number which the currently opened file is attached to,
\var{iscan} is the current scan, \var{ispecerr} indicates something is going
wrong.
\var{ncolumns} is the number of columns of data.
The parameters \param{LINELENGTH} and \param{WORDLENGTH} specify the length
of some character strings and
\param{NWORDS} and \param{NLINES} dimension arrays for us.
\param{LINELENGTH} is the maximum length of a line (distance between
carriage returns) which we allow.
\var{specsfilename} will hold the name that SPEC was using for the name
of the file and \var{specsdate} will hold the date the file
inwas created on (when the experiment was originally carried out).
This information can be inserted into datafile titles if nothing else is
provided.
We think (but are not sure) that longer lines will be truncated if read into
a string which is too short.
The integer array \var{headerwords} makes a note of how many
words are found on the
\#O lines and \param{NLINES} is the maximum number of \#O lines we can handle.
Those words will be stored in the array \var{words}.
\var{wordvalues} will hold the numbers corresponding to the words held in
\var{words}.
\var{filnam} is the name of the SPEC file and line is the
contents of the most recently read line.
The scan type (turboscan, ascan, whatever) will be stored in \var{scantype}.
Finally the column labels (from the \#L) line are stored in \var{columnlabels}.
\subsection{File opening}
A file is needed before anything can be done with it so we start by trying
to get hold of one.
At this point we'll need to know the filename which the user is going to have
to specify if it wasn't already supplied.
Calling the routine when a file is already open causes that file to be
closed and the current one to be opened (so only one file will ever be open
at once from here).
Normally the filename will be supplied from the command line, and the
main program should ensure it is already filled in for us.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap2}\raggedright\small
\NWtarget{nuweb4}{} $\langle\,${\it getfile}\nobreak\ {\footnotesize {4}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine getfile@\\
\mbox{}\verb@ logical :: od@\\
\mbox{}\verb@! See if there is already a file open, if there is then close it@\\
\mbox{}\verb@ inquire(unit=iunit,opened=od)@\\
\mbox{}\verb@ if(od)close(iunit)@\\
\mbox{}\verb@! filnam must already be filled in, or we'll try to open garbage@\\
\mbox{}\verb@1 open(unit=iunit,file=filnam,form='formatted', &@\\
\mbox{}\verb@ & access='sequential', status='old',iostat=ispecerr)@\\
\mbox{}\verb@ if(ispecerr.ne.0) then@\\
\mbox{}\verb@ write(*,'(3a,i6)')'Problem with file ',filnam(1:len_trim(filnam)),&@\\
\mbox{}\verb@ & ' iostat=',ispecerr@\\
\mbox{}\verb@ write(*,'(a)')'Please supply a valid SPEC file name'@\\
\mbox{}\verb@ read(*,lnfmt)filnam@\\
\mbox{}\verb@ goto 1@\\
\mbox{}\verb@ endif @\\
\mbox{}\verb@ topofscan=.false.@\\
\mbox{}\verb@! Supply defaults in case dumb specfile has no #L@\\
\mbox{}\verb@ ncolumns=14 ; columnlabels(1)='2_theta'; columnlabels(2)='Epoch' @\\
\mbox{}\verb@ columnlabels(3)='Seconds'; columnlabels(4)='MA0' @\\
\mbox{}\verb@ columnlabels(5)='MA1'; columnlabels(6)='MA2' @\\
\mbox{}\verb@ columnlabels(7)='MA3'; columnlabels(8)='MA4' @\\
\mbox{}\verb@ columnlabels(9)='MA5'; columnlabels(10)='MA6' @\\
\mbox{}\verb@ columnlabels(11)='MA7'; columnlabels(12)='MA8' @\\
\mbox{}\verb@ columnlabels(13)='Monitor'; columnlabels(14)='Fluo det'@\\
\mbox{}\verb@ i0=ichar('0');i9=ichar('9')@\\
\mbox{}\verb@ id=ichar('.');isp=ichar(' ');im=ichar('-')@\\
\mbox{}\verb@ return@\\
\mbox{}\verb@ end subroutine getfile @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
Default values are provided for the number of data columns, \var{ncolumns},
and their labels, in \var{columnlabels}. The values \var{i0}, \var{i9},
\var{id}, \var{is} and \var{im} are used in the \code{rdnums}
routine, and initialised here as we expect this will
always be called before we try to read any numbers.
This will need to be modified or replaced if it sits behind a gui, or perhaps
it should just send a unit number to \mod{useroptions}?
\subsection{Finding a specific scan}
Having gotten a file opened we'll want to read it.
Usually we'll go after a particular scan so here's a routine for doing
just that, provided the scans are in order in the SPEC file.
It positions the file at the top of the data beginning of the scan.
Any header information will already be read in by the \code{readheader}
routine, so there is no need to call \code{readheader} as well as
\code{findscan}. (That would actually cause problems).
There is an implicit assumption that scans will always be in the file
in ascending order and that we will always access them in ascending
order. If this is not the case we can add some kind of SPEC file preprocessor
which puts the scans in order, which is more suited to some kind of script.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap3}\raggedright\small
\NWtarget{nuweb5}{} $\langle\,${\it findscan}\nobreak\ {\footnotesize {5}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine findscan(n)@\\
\mbox{}\verb@ integer(kind=4),intent(in)::n@\\
\mbox{}\verb@ character(len=7) :: rd ! enough space for yes, no and unknown@\\
\mbox{}\verb@ inquire(unit=iunit,read=rd) ! Can we read the file?@\\
\mbox{}\verb@ if(rd(1:3).ne.'YES') call getfile ! If not go open it @\\
\mbox{}\verb@ if(iscan.eq.n .and. topofscan) return@\\
\mbox{}\verb@ if(iscan.gt.n) rewind(n) ! @\\
\mbox{}\verb@2 call readheader ! in case top of file @\\
\mbox{}\verb@ if(iscan.ne.n)then ! work through file to find #@\\
\mbox{}\verb@1 read(iunit,lnfmt,err=10,end=10)line@\\
\mbox{}\verb@ if(line(1:1).eq.'#') goto 2@\\
\mbox{}\verb@ go to 1@\\
\mbox{}\verb@ else@\\
\mbox{}\verb@! Scan found - check scantype elsewhere - not a job for specfile module@\\
\mbox{}\verb@ return@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@10 write(*,*)@\\
\mbox{}\verb@ write(*,'(a,i5,a)')'Scan ',n,' not found in file '// &@\\
\mbox{}\verb@ & filnam(1:len_trim(filnam))@\\
\mbox{}\verb@ ispecerr=-2@\\
\mbox{}\verb@ return@\\
\mbox{}\verb@ end subroutine findscan @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
The routine returns sets \var{ispecerr} to -2 if the scan is not
found in the file.
This is to distinguish from the end of a scan, which is currently
caught by setting \var{ispecerr} to -1.
Note that the routine uses the \code{readheader} routine to interpret
any header lines found, rather than trying to read them itself.
The logical \var{topofscan} is set by \code{readheader} when it gets to
a data line after a set of header lines.
These complications are just to prevent us ever needing to read the
same line twice.
\subsection{Reading and interpreting file headers}
It will be necessary to read and intelligently interpret the various bits
of header information which go with each scan.
At the very least we need to check the kind of scan and the header information
for each of the columns.
If the detector offsets ever get stored in the specfile this is a
good time to dump them on a temp.res file or something.
Probably strain scanning experiments will want various motor positions
recorded somewhere, so we just pull everything possible out here.
For interpreting the header we need some information about the meaning of the
various \# labels in the header. This has been alluded to earlier, but here is
a definitive list.
\begin{itemize}
\item \#F filename
\item \#E A number - not sure what it is!
\item \#D The date
\item \#C The creator of the file
\item \#On Index of the labels which will correspond to following \#P cards
\item \#Sn Scan number and information
\item \#G Diffractometer geometry for spec fourc mode (see link below)
\item \#Q Four circle variables. h,k,l = q[0,1,2], wavelength=q[3]
\item \#Pn values corresponding to motor names in \#O
\item \#N number of columns in output data
\item \#L names of columns in output data (should be \#N columns)
\item blank, the data themselves
\item blank line - end of a scan?
\item scan ended by control-c ? Did I imagine this ?
\end{itemize}
For the rest of \#Q see: http://www.certif.com/spec\_manual/fourc\_4\_9.html
The \code{readheader} routine interprets these lines.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap4}\raggedright\small
\NWtarget{nuweb6}{} $\langle\,${\it readheader}\nobreak\ {\footnotesize {6}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine readheader@\\
\mbox{}\verb@ character(len=1) :: letter@\\
\mbox{}\verb@ character(len=128) :: motor@\\
\mbox{}\verb@ integer(kind=4) :: i,j@\\
\mbox{}\verb@ real(kind=8) :: a@\\
\mbox{}\verb@ integer(kind=4),parameter :: four=4@\\
\mbox{}\verb@ if(line(1:1).eq.'#') goto 2 ! if already on header, don't skip@\\
\mbox{}\verb@1 read(iunit,lnfmt,end=100)line@\\
\mbox{}\verb@2 if(line(1:1).eq.'#') then; letter=line(2:2); topofscan=.true.@\\
\mbox{}\verb@ select case((letter))@\\
\mbox{}\verb@ case('O')@\\
\mbox{}\verb@ read(line(3:3),'(i1)')i@\\
\mbox{}\verb@ call split(line(4:LINELENGTH),words(:,i+1),WORDLENGTH,NWORDS,j)@\\
\mbox{}\verb@ headerwords(i+1)=j@\\
\mbox{}\verb@ case('P')@\\
\mbox{}\verb@ read(line(3:3),'(i1)')i@\\
\mbox{}\verb@ call rdnums(four,NWORDS,wordvalues(:,i+1))@\\
\mbox{}\verb@ case('Q')@\\
\mbox{}\verb@! Copy the Q from the header into our specfile module@\\
\mbox{}\verb@ read(line(3:LINELENGTH),*,end=1, err=1) Q@\\
\mbox{}\verb@ case('N')@\\
\mbox{}\verb@ read(line(3:len_trim(line)),*,end=1)ncolumns@\\
\mbox{}\verb@ case('L') ! signals end of header, bug out here !@\\
\mbox{}\verb@ call split(line(3:LINELENGTH),columnlabels,WORDLENGTH,NWORDS,i)@\\
\mbox{}\verb@ if(ncolumns.ne.i) &@\\
\mbox{}\verb@ & write(*,'(a,i5)')'error reading header for scan',iscan @\\
\mbox{}\verb@ case('S')@\\
\mbox{}\verb@ read(line(3:len_trim(line)),*,end=1)iscan,scantype,motor, &@\\
\mbox{}\verb@ & scanstart, scanend@\\
\mbox{}\verb@ case('F')@\\
\mbox{}\verb@ specsfilename=line(3:len_trim(line)) ! get original filename@\\
\mbox{}\verb@ case('D')@\\
\mbox{}\verb@ read(line(3:len_trim(line)),'(a256)',end=1)specsdate@\\
\mbox{}\verb@ case('C'); continue ! comments@\\
\mbox{}\verb@ case('G'); continue ! no idea - always zero@\\
\mbox{}\verb@ case('E'); continue ! epoch - we don't care what it was for now.@\\
\mbox{}\verb@ case default@\\
\mbox{}\verb@ end select@\\
\mbox{}\verb@ else ! Not a # line, so assume end of header@\\
\mbox{}\verb@ read(line,*,err=1,end=1)a ! catch blank lines in header@\\
\mbox{}\verb@ topofscan=.true. !!! Must be able to read a number !@\\
\mbox{}\verb@ return !!! escapes from routine here !!!!!@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ goto 1@\\
\mbox{}\verb@100 ispecerr=-1; return; end subroutine readheader @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
A short routine to split a line into a series of words was needed in the
\code{readheader} routine and might generally come in useful, so here it is.
Note that the gap between words has to be at least two spaces to allow motor
names which contain a space.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap5}\raggedright\small
\NWtarget{nuweb7a}{} $\langle\,${\it split}\nobreak\ {\footnotesize {7a}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine split(instring,outstrings,lenout,n,i)@\\
\mbox{}\verb@ integer(kind=4),intent(in) :: lenout,n@\\
\mbox{}\verb@ character(len=*),intent(in) :: instring@\\
\mbox{}\verb@ character(len=lenout),dimension(n),intent(out) :: outstrings@\\
\mbox{}\verb@ integer(kind=4),intent(out) :: i@\\
\mbox{}\verb@ integer(kind=4) :: j,k,l@\\
\mbox{}\verb@ j=1; k=1@\\
\mbox{}\verb@ do i=1,len_trim(instring) ! hope len > len_trim or array overstep@\\
\mbox{}\verb@ if(instring(i:i+1).ne.' ')then@\\
\mbox{}\verb@ outstrings(j)(k:k)=instring(i:i)@\\
\mbox{}\verb@ if(k.ne.1 .and. k.lt.lenout) k=k+1@\\
\mbox{}\verb@ if(instring(i:i).ne.' '.and.k.eq.1)k=k+1@\\
\mbox{}\verb@ else @\\
\mbox{}\verb@ if(k.gt.1)then ; do l=k,lenout@\\
\mbox{}\verb@ outstrings(j)(l:l)=' ' ! blank pad end of string@\\
\mbox{}\verb@ enddo ; j=j+1; k=1; if(j.gt.n)exit ; endif@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ enddo @\\
\mbox{}\verb@! Blank pad last string if necessary@\\
\mbox{}\verb@ if(k.gt.1)then@\\
\mbox{}\verb@ do l=k,lenout;outstrings(j)(l:l)=' ';enddo@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ i=j; return; end subroutine split @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
A quite involded routine which was needed in \code{readheader}, to convert a
string containing a space separated list of numbers into an array.
This same routine is used from the \code{getdata} routine.
It is designed to be significantly faster than a simplistic
fortran \code{read(line,*) vars}, at least on some platforms.
Briefly, it parses along the string, ignoring any leading white space and then
inteprets the characters it finds into and array of numbers.
The only characters a line is allowed to contain (currently) are therefore
1,2,3,4,5,6,7,8,9,0,\verb*+ +,. and -.
Anything else causes \var{ispecerr} to be set to -1, which normally signals the
end of the file (FIXME - add E formats?).
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap6}\raggedright\small
\NWtarget{nuweb7b}{} $\langle\,${\it rdnums}\nobreak\ {\footnotesize {7b}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine rdnums(ic,n,values)@\\
\mbox{}\verb@! Placed here in a subroutine in case of formatting or error handling problems@\\
\mbox{}\verb@ integer(kind=4),intent(in) :: n, ic@\\
\mbox{}\verb@ real(kind=8),dimension(n),intent(inout) :: values@\\
\mbox{}\verb@ if(line(ic:ic).eq.'#')then; ispecerr=-1; return; endif@\\
\mbox{}\verb@ read(line(ic:len_trim(line)),*,err=10,end=20)values(1:ncolumns)@\\
\mbox{}\verb@ goto 100@\\
\mbox{}\verb@10 write(*,*)'Error reading line, looking for ',ncolumns,' values'@\\
\mbox{}\verb@ write(*,*)line(ic:len_trim(line))@\\
\mbox{}\verb@ write(*,*)values@\\
\mbox{}\verb@20 continue@\\
\mbox{}\verb@100 return@\\
\mbox{}\verb@ end subroutine rdnums @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap7}\raggedright\small
\NWtarget{nuweb8}{} $\langle\,${\it junk1}\nobreak\ {\footnotesize {8}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@@\\
\mbox{}\verb@ integer(kind=4) :: i, j, k, l; real(kind=8) :: t@\\
\mbox{}\verb@ integer(kind=4) :: ii0,ii9,iid,iisp,iim@\\
\mbox{}\verb@ logical :: pastdot, neg , innum@\\
\mbox{}\verb@ k=1; t=0.1d0; l=0 ; values=0.0d0@\\
\mbox{}\verb@ ii0=ichar('0');ii9=ichar('9');iisp=ichar(' ');iid=ichar('.');iim=ichar('-')@\\
\mbox{}\verb@ innum=.false. ; pastdot=.false. ; neg=.false.@\\
\mbox{}\verb@ do i=ic,LINELENGTH ! strip leading spaces@\\
\mbox{}\verb@ if(ichar(line(i:i)).ne.iisp)then; l=i; exit ; endif; enddo@\\
\mbox{}\verb@ if(l.eq.0)then; ispecerr=-1 ; return ; endif@\\
\mbox{}\verb@ do i=l,LINELENGTH ! instring is always line@\\
\mbox{}\verb@ j=ichar(line(i:i))@\\
\mbox{}\verb@ if(j.le.ii9 .and. j.gt.ii0) then ! digit@\\
\mbox{}\verb@ innum=.true.@\\
\mbox{}\verb@ if(pastdot)then@\\
\mbox{}\verb@ values(k)=values(k) + real(j-ii0,8)*t ; t=t/10.0d0@\\
\mbox{}\verb@ else@\\
\mbox{}\verb@ values(k)=values(k)*10.0d0 + real(j-ii0,8)@\\
\mbox{}\verb@ endif; cycle ;endif !!! next char !!!@\\
\mbox{}\verb@ if(j.eq.iisp)then ! space@\\
\mbox{}\verb@ if(.not.innum)cycle !!! next char !!!@\\
\mbox{}\verb@ if(neg)values(k)=-values(k)@\\
\mbox{}\verb@ if(k.lt.n) then@\\
\mbox{}\verb@ k=k+1@\\
\mbox{}\verb@ else ; return; @\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ t=0.1d0; pastdot=.false.; neg=.false.; innum=.false.@\\
\mbox{}\verb@ cycle; @\\
\mbox{}\verb@ endif !!! next char !!!@\\
\mbox{}\verb@ if(j.eq.ii0)then ; innum=.true.@\\
\mbox{}\verb@ if(pastdot) then ; t=t/10.0d0 ; else@\\
\mbox{}\verb@ values(k)=values(k)*10.0d0 ; endif ; cycle; endif @\\
\mbox{}\verb@ if(j.eq.iid)then; innum=.true.; pastdot=.TRUE. ; cycle ; endif@\\
\mbox{}\verb@ if(j.eq.iim)then; innum=.true.; neg=.TRUE. ; cycle ; endif@\\
\mbox{}\verb@ if(line(ic:ic).eq.'#')then; ispecerr=-1; return; endif@\\
\mbox{}\verb@ if(line(i:i).eq.'e' .or. line(i:i).eq.'E')then @\\
\mbox{}\verb@ values=0.0d0 ! give up and use fortran read !@\\
\mbox{}\verb@! write(*,*)'About to call fortran read'@\\
\mbox{}\verb@ read(line(ic:len_trim(line)),err=10,end=20)values(1:ncolumns)@\\
\mbox{}\verb@ goto 100 ! bug out after successful 'e' format read@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@! error reading line@\\
\mbox{}\verb@ write(*,*)'Could not interpret character... line was:'@\\
\mbox{}\verb@ write(*,'(a)')line(1:len_trim(line))@\\
\mbox{}\verb@ write(*,'(a,i4,a,i4)')line(ic:ic),j,'at position',ic@\\
\mbox{}\verb@! write(*,*)'i0',ii0,'ichar(0)=',ichar('0'),'i9=',ii9,'=',ichar('9')@\\
\mbox{}\verb@10 continue@\\
\mbox{}\verb@ enddo @\\
\mbox{}\verb@20 continue@\\
\mbox{}\verb@100 return; end subroutine rdnums @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item {\NWtxtMacroNoRef}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap8}\raggedright\small
\NWtarget{nuweb9a}{} $\langle\,${\it junk}\nobreak\ {\footnotesize {9a}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@@\\
\mbox{}\verb@10 write(*,'(a)')'rdnums couldn''t interpret the line'@\\
\mbox{}\verb@ write(*,*)'ic=',ic,'n=',n,'j=',j@\\
\mbox{}\verb@ write(*,'(a)')line(ic:len_trim(line))@\\
\mbox{}\verb@ values=0.0d0@\\
\mbox{}\verb@ ispecerr=-1; @\\
\mbox{}\verb@ return ! Give error message if bad character@\\
\mbox{}\verb@! end of line reached - OK as we were able to get line from file@\\
\mbox{}\verb@20 return ! with hopefully first few values filled in@\\
\mbox{}\verb@ enddo; return; end subroutine rdnums @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item {\NWtxtMacroNoRef}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
This horribly complicated routine might later get replaced with something
else. It gave a factor of more than 2 (55 seconds to 18 seconds) for
reading a 50MB file.
To make head or tail of what it is doing you would probably want to expand
all the semicolons to put the code on separate lines.
It has been tested and it thought to be correct however - it's fairly
straightforward to watch under a debugger to see if it reads a string
correctly. Note that \var{line} should already be filled in, this routine
reads whatever is in that variable.
Obtaining the information from the file headers, indexed by keyword requires
a small function which sifts through the headers which have been read in and
then gets the corresponding value.
This is to be called after the header is read in, if required.
It will be able to supply the initial 2$\theta$ position and the positions of
any motors of interest (like X trans, Y trans, Chi etc).
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap9}\raggedright\small
\NWtarget{nuweb9b}{} $\langle\,${\it getheadervalue}\nobreak\ {\footnotesize {9b}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ real(kind=8) function getheadervalue(string)@\\
\mbox{}\verb@ character(len=*),intent(in) :: string@\\
\mbox{}\verb@ integer(kind=4) :: i, j, k@\\
\mbox{}\verb@ k=len(string); getheadervalue=0.0@\\
\mbox{}\verb@ do i=1,NWORDS; do j=1,NLINES@\\
\mbox{}\verb@ if(string(1:k).eq.words(i,j)(1:k))then@\\
\mbox{}\verb@ getheadervalue=wordvalues(i,j); return @\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ enddo; enddo@\\
\mbox{}\verb@ return; end function getheadervalue @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
Usage is then simply \code{variable = getheadervalue('X translation')}.
We might add a logfile option to binem to dump out selected motors - will
need to do something quite general about spaces in motor names here, they
cannot be supplied as command line arguments without quoting.
Understanding the contents of the data columns when they are read in means we
need to be able to ask for the column label of a particular counter.
The next routine goes through the \var{columnlabels} array which was filled
in during \code{readheader} to work out which column a particular label
corresponds to.
It returns -1 if the column is not there (use for example to work out if we
have recorded a particular counter value).
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap10}\raggedright\small
\NWtarget{nuweb10a}{} $\langle\,${\it whichcolumn}\nobreak\ {\footnotesize {10a}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ integer(kind=4) function whichcolumn(string)@\\
\mbox{}\verb@! Interprets the #L line information@\\
\mbox{}\verb@ character(len=*),intent(in) :: string@\\
\mbox{}\verb@ integer(kind=4) :: i, j@\\
\mbox{}\verb@ j=len(string) ; whichcolumn=-1@\\
\mbox{}\verb@ do i=1,NWORDS@\\
\mbox{}\verb@ if(string(1:j).eq.columnlabels(i)(1:j))then@\\
\mbox{}\verb@ whichcolumn=i ; return@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ enddo@\\
\mbox{}\verb@ return; end function whichcolumn @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
\subsection{Reading lines of data}
At last we come to the most interesting thing, which is to read in a line of
data.
Arguments to a subroutine are used to pass the information back.
The calling routine is going to have to specify how many numbers to read and
supply an array to put them into
All of the glitch elimination, zinger elimination, detection of closed
shutters and whatnot is going to happen elsewhere - here we only send back
some numbers and set \var{ispecerr} to a non zero value if something went
wrong.
The caller should check this value and any data returned from here.
A futuristic situation would be to have a version of this routine
which can pass back the data as it arrives from the instrumentation, along
some kind of unix pipe thing.
A quick test on the amber cluster (linux computers at ESRF) showed that
if a file is appended with \code{cat >> file} then a fortran program does
indeed see the extra lines even if they appear after it starts running.
However, if the top of the file is modified (so the fortran program would lose
its place) then it seems to work off of a copy of the file in the state
it was in when it was opened. These details can be worked out as and when
we try for an online binning program.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap11}\raggedright\small
\NWtarget{nuweb10b}{} $\langle\,${\it getdata}\nobreak\ {\footnotesize {10b}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ subroutine getdata(a,n)@\\
\mbox{}\verb@! A is an array, dimensioned by the number of data items to be expected@\\
\mbox{}\verb@ integer(kind=4),intent(in) :: n@\\
\mbox{}\verb@ real(kind=8),dimension(n),intent(out) :: a@\\
\mbox{}\verb@ integer(kind=4), parameter :: one = 1@\\
\mbox{}\verb@ if(topofscan)then; topofscan=.false. ; else@\\
\mbox{}\verb@ read(iunit,lnfmt,err=10,end=10)line ! goes via line@\\
\mbox{}\verb@ topofscan=.false.@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ call rdnums(one,n,a) ! reads from line@\\
\mbox{}\verb@ if(line(1:1).eq.'#')then; ispecerr=-1 ;goto 100 ; endif@\\
\mbox{}\verb@ go to 100 @\\
\mbox{}\verb@10 ispecerr=-1 ! End of file@\\
\mbox{}\verb@100 return ; end subroutine getdata @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11a}{11a}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
The \code{getdata} routine is just calling the \code{rdnums} routine to
interpret the line of numbers.
This assume all numbers can sensibly be read into a floating point format (so
integers are converted to reals).
\subsection{A module for all of that specfile stuff}
Just placing all the variables and subroutines together in a module for easy
access elsewhere. Apart from trying out a fortran 90 feature the author was
curious about, this is supposed to be an aid to structured programming.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap12}\raggedright\small
\NWtarget{nuweb11a}{} $\langle\,${\it specfiles}\nobreak\ {\footnotesize {11a}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ module specfiles@\\
\mbox{}\verb@@\hbox{$\langle\,${\it SPECVARS}\nobreak\ {\footnotesize \NWlink{nuweb3}{3}}$\,\rangle$}\verb@@\\
\mbox{}\verb@ contains@\\
\mbox{}\verb@@\hbox{$\langle\,${\it getfile}\nobreak\ {\footnotesize \NWlink{nuweb4}{4}}$\,\rangle$}\verb@ @\\
\mbox{}\verb@@\hbox{$\langle\,${\it findscan}\nobreak\ {\footnotesize \NWlink{nuweb5}{5}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it getdata}\nobreak\ {\footnotesize \NWlink{nuweb10b}{10b}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it readheader}\nobreak\ {\footnotesize \NWlink{nuweb6}{6}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it getheadervalue}\nobreak\ {\footnotesize \NWlink{nuweb9b}{9b}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it whichcolumn}\nobreak\ {\footnotesize \NWlink{nuweb10a}{10a}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it rdnums}\nobreak\ {\footnotesize \NWlink{nuweb7b}{7b}}$\,\rangle$}\verb@@\\
\mbox{}\verb@@\hbox{$\langle\,${\it split}\nobreak\ {\footnotesize \NWlink{nuweb7a}{7a}}$\,\rangle$}\verb@@\\
\mbox{}\verb@ end module specfiles @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item \NWtxtMacroRefIn\ \NWlink{nuweb11b}{11b}\NWlink{nuweb12}{, 12}\NWlink{nuweb13}{, 13}\NWlink{nuweb14b}{, 14b}\NWlink{nuweb15}{, 15}\NWlink{nuweb16}{, 16}\NWlink{nuweb31a}{, 31a}\NWlink{nuweb32}{, 32}\NWlink{nuweb40a}{, 40a}\NWlink{nuweb57b}{, 57b}\NWlink{nuweb86b}{, 86b}\NWlink{nuweb89}{, 89}\NWlink{nuweb91}{, 91}.
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
So in theory any programmer can come along in the future and if they add the
line \code{use specfiles} to their programs then they'll have access to all
those goodies. The file \code{specfiles.f90} should be compilable into
a library (type thing).
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap13}\raggedright\small
\NWtarget{nuweb11b}{} \verb@"specfiles.f90"@\nobreak\ {\footnotesize {11b}}$\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@@\hbox{$\langle\,${\it specfiles}\nobreak\ {\footnotesize \NWlink{nuweb11a}{11a}}$\,\rangle$}\verb@ @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
\subsection{Test program}
In order to check the routines in this section are working we have a simple
program which gets a file name from the command line (or prompts)
and dumps out header and column information on stdout.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap14}\raggedright\small
\NWtarget{nuweb12}{} \verb@"testspec.f90"@\nobreak\ {\footnotesize {12}}$\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@@\hbox{$\langle\,${\it specfiles}\nobreak\ {\footnotesize \NWlink{nuweb11a}{11a}}$\,\rangle$}\verb@@\\
\mbox{}\verb@ program testspec@\\
\mbox{}\verb@ use specfiles@\\
\mbox{}\verb@ integer(kind=4) :: i,j@\\
\mbox{}\verb@ real(kind=8),allocatable :: data(:), prevdata(:)@\\
\mbox{}\verb@ real :: time1, time2@\\
\mbox{}\verb@ call getarg(1,filnam)@\\
\mbox{}\verb@ call getfile@\\
\mbox{}\verb@ call cpu_time(time1)@\\
\mbox{}\verb@ call readheader ! Get's the header at the top of the file@\\
\mbox{}\verb@2 write(*,'(a5,i5)')'Scan ',iscan@\\
\mbox{}\verb@ write(*,'(a)')'Header words and values'@\\
\mbox{}\verb@ do j=1,NLINES @\\
\mbox{}\verb@ if(headerwords(j).gt.0)then@\\
\mbox{}\verb@ do i=1,headerwords(j) @\\
\mbox{}\verb@ write(*,'(a20,f16.8)')words(i,j),wordvalues(i,j)@\\
\mbox{}\verb@ enddo@\\
\mbox{}\verb@ endif ! headerwords(j.gt.0)@\\
\mbox{}\verb@ enddo@\\
\mbox{}\verb@ write(*,'(a,i5)')'Number of columns = ', ncolumns@\\
\mbox{}\verb@ if(.not.allocated(data))allocate(data(ncolumns))@\\
\mbox{}\verb@ if(.not.allocated(prevdata))allocate(prevdata(ncolumns))@\\
\mbox{}\verb@ if(size(data).ne.ncolumns)then@\\
\mbox{}\verb@ deallocate(data); deallocate(prevdata)@\\
\mbox{}\verb@ allocate(data(ncolumns));allocate(prevdata(ncolumns))@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ j=0@\\
\mbox{}\verb@1 prevdata=data ! Start of loop through reading data@\\
\mbox{}\verb@ call getdata(data,ncolumns)@\\
\mbox{}\verb@ j=j+1@\\
\mbox{}\verb@ if(ispecerr.eq.0)goto 1 ! loop to end of scan@\\
\mbox{}\verb@ ispecerr=0 ! OK, we reached the end of the scan@\\
\mbox{}\verb@ write(*,'(a)')'Last data line in scan'@\\
\mbox{}\verb@ do i=1,ncolumns@\\
\mbox{}\verb@ write(*,'(a20,f16.8)')columnlabels(i),prevdata(i)@\\
\mbox{}\verb@ enddo@\\
\mbox{}\verb@ write(*,'(a,i10)')'Number of data points = ',j@\\
\mbox{}\verb@ call findscan(iscan+1) ! get the next scan (assumes sequential)@\\
\mbox{}\verb@ if(ispecerr.eq.0)goto 2 ! process next scan@\\
\mbox{}\verb@ deallocate(data); deallocate(prevdata)@\\
\mbox{}\verb@ call cpu_time(time2)@\\
\mbox{}\verb@ write(*,*)@\\
\mbox{}\verb@ write(*,'(a,f8.4,a)')'Time to read file = ',time2-time1,' seconds'@\\
\mbox{}\verb@ write(*,'(a)')'Thats all folks!'@\\
\mbox{}\verb@ end program testspec @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
This program could form the basis of a rebinning program - all it needs to do
as extra work is to determine the user options for which scan to bin, which
stepsize etc and then pass the data to the rebinning, summing and output
routines.
It could also be modified to write out summary information about the scan
(total counts) or to pull out information about temperatures and so on.
During testing it was noted that there is a factor 5 in the speed that the
program runs if the line of data is actually read in, or just checked
to see if it the second character on the line is a digit.
This suggests we might try to optimise the rdnums subroutine.
(FIXME) Some more work will be necessary (somewhere) to catch corrupted
specfiles, and allow files to be read where scans are not in order.
\subsection{Utility programs}
Some useful programs for extracting data from SPEC files are given here.
Plotit just pulls out a column of data from the SPEC file and creates
an mtv format file for display using the plotmtv program.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap15}\raggedright\small
\NWtarget{nuweb13}{} \verb@"plotit.f90"@\nobreak\ {\footnotesize {13}}$\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@@\hbox{$\langle\,${\it specfiles}\nobreak\ {\footnotesize \NWlink{nuweb11a}{11a}}$\,\rangle$}\verb@@\\
\mbox{}\verb@ program plotit@\\
\mbox{}\verb@ use specfiles@\\
\mbox{}\verb@ integer(kind=4) :: j@\\
\mbox{}\verb@ integer(kind=4) :: nscan, ncol@\\
\mbox{}\verb@ character(LEN=WORDLENGTH) :: label@\\
\mbox{}\verb@ real(kind=8),allocatable :: data(:), prevdata(:)@\\
\mbox{}\verb@ real :: time1, time2@\\
\mbox{}\verb@@\hbox{$\langle\,${\it cmdline}\nobreak\ {\footnotesize \NWlink{nuweb14a}{14a}}$\,\rangle$}\verb@@\\
\mbox{}\verb@ call cpu_time(time1)@\\
\mbox{}\verb@ if(.not.allocated(data))allocate(data(ncolumns))@\\
\mbox{}\verb@ if(.not.allocated(prevdata))allocate(prevdata(ncolumns))@\\
\mbox{}\verb@ if(size(data).ne.ncolumns)then@\\
\mbox{}\verb@ deallocate(data); deallocate(prevdata)@\\
\mbox{}\verb@ allocate(data(ncolumns));allocate(prevdata(ncolumns))@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ j=0@\\
\mbox{}\verb@ write(*,'(a)')'$ DATA = CURVE2D'@\\
\mbox{}\verb@ write(*,'(3a)')'% xlabel = "',columnlabels(1),'"'@\\
\mbox{}\verb@ write(*,'(3a)')'% ylabel = "',columnlabels(ncol),'"'@\\
\mbox{}\verb@ write(*,'(a)')'% linetype = 1 markertype = 2'@\\
\mbox{}\verb@ write(*,'(a)')'% linecolor = 3 markercolor = 4'@\\
\mbox{}\verb@1 prevdata=data ! Start of loop through reading data@\\
\mbox{}\verb@ call getdata(data,ncolumns)@\\
\mbox{}\verb@ j=j+1@\\
\mbox{}\verb@ if(ispecerr.eq.0)then@\\
\mbox{}\verb@ write(*,*)data(1),data(ncol)@\\
\mbox{}\verb@ goto 1@\\
\mbox{}\verb@ endif@\\
\mbox{}\verb@ deallocate(data); deallocate(prevdata)@\\
\mbox{}\verb@ call cpu_time(time2)@\\
\mbox{}\verb@ write(*,*)'$ END'@\\
\mbox{}\verb@! write(*,'(a,f8.4,a)')'Time to read file = ',time2-time1,' seconds'@\\
\mbox{}\verb@! write(*,'(a)')'Thats all folks!'@\\
\mbox{}\verb@100 end program plotit @{\NWsep}
\end{list}
\vspace{-1.5ex}
\footnotesize
\begin{list}{}{\setlength{\itemsep}{-\parsep}\setlength{\itemindent}{-\leftmargin}}
\item{}
\end{list}
\end{minipage}\vspace{4ex}
\end{flushleft}
The command line stuff is generally the same from some of the other programs
below, so we separate it off here.
\begin{flushleft} \small
\begin{minipage}{\linewidth}\label{scrap16}\raggedright\small
\NWtarget{nuweb14a}{} $\langle\,${\it cmdline}\nobreak\ {\footnotesize {14a}}$\,\rangle\equiv$
\vspace{-1ex}
\begin{list}{}{} \item
\mbox{}\verb@ call getarg(1,filnam)@\\
\mbox{}\verb@ call getfile@\\
\mbox{}\verb@ call getarg(2,line)@\\
\mbox{}\verb@ read(line,*,err=10,end=10)nscan@\\
\mbox{}\verb@ call findscan(nscan)@\\
\mbox{}\verb@ call getarg(3,label)@\\
\mbox{}\verb@ read(label,*,err=11,end=11)ncol@\\
\mbox{}\verb@ goto 20@\\
\mbox{}\verb@! allow labels as well as numbers for column headings@\\