Ticket #562: asym_rgxx_nexus.for

File asym_rgxx_nexus.for, 13.2 KB (added by Nick Draper, 11 years ago)
Line 
1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2C
3      PROGRAM ASYM_RG
4C     ---------------
5C
6C  9-FEB-97 C A Scott   Corrected typo in error calc in INTSYM
7C                       Effect of typo was to underestimate
8C                       error by a factor of 1.4
9C
10C 18-MAY-1998 J S Lord  Version for 16 Red & 16 Green histograms on EMU
11C
12C    NOV-1998 S P Cottrell  Version for 32 Red & 32 Green histograms on muSR
13C
14C  7-JAN-1999 S P Cottrell  Added run numbers in titles on graph, fixed size
15C                           of graph components and added facility to
16C                           create a PGPLOT.PS file immediately after
17C                           viewing plot
18C
19C 29-APR-2006 S P Cottrell Reads only NeXus files
20C
21C 4-OCT-2006  S P Cottrell Can dig scans out of multiperiod data files
22C
23C
24      REAL      X(200),YI(200),YD(200),C(256000)
25      REAL      F(2,2048,200),B(2,2048,200)
26      REAL      EI(200),ED(200),HP,H
27      LOGICAL   AGAIN,LMORE
28      INTEGER   LRG,RP,GP
29      CHARACTER TITLE*60
30        character*4 beamstr
31        integer beam
32C
33      write(*,*)'Multiperiod scans'
3477      write(*,78)
3578      format(' Beamline? ',$)
36        read(*,79)beamstr
3779      format(A4)
38        if(beamstr(1:3).eq.'EMU' .or. beamstr(1:3).eq.'emu') then
39          beam=1
40        else if(beamstr(1:4).eq.'MUSR' .or. beamstr(1:4).eq.'musr') then
41          beam=2
42        else if(beamstr(1:4).eq.'DEVA' .or. beamstr(1:4).eq.'deva') then
43          beam=3
44        else
45          goto 77
46        endif
47      WRITE(*,*)
48
49C      CALL ASKSIZ(SIZE,IWIDTH)
50C     We probably don't need this facility, so fix character size and
51C     linewidth
52      SIZE = 1.25
53      IWIDTH = 1
54
55      WRITE(*,*)
56      K=0
57   1  WRITE(*,*)
58      WRITE(*,100)
59 100  FORMAT(' >  First and last run number ?  : ',$)
60      READ(*,*,ERR=1) ISTART,IEND
61      IF (ISTART.LT.1 .OR. IEND.LT.ISTART) GOTO 1
62      NRUNS=IEND-ISTART+1
63      IF ((K+NRUNS).GT.200) THEN
64        WRITE(*,*)' No more than 200 runs allowed!'
65        GOTO 1
66      ENDIF
67      DO 10 IRUN=ISTART,IEND
68        K=K+1
69        IF (IRUN.EQ.ISTART) THEN
70           CALL GETFIL(IRUN,NT,NH,H,C,DXTIME,NFRAMS,
71     +          TITLE,LRG,.TRUE.,beam)
72           PRINT *,' '
73           PRINT *,' Number of periods: ', LRG
74           PRINT *,' Total number of histograms: ', NH
75           HP = NH/LRG
76           PRINT *,' Histograms in each period: ', HP
77
78  150      PRINT *,' Enter Period for ''Red'' data'
79           READ (*,*) RP
80           PRINT *,' Enter Period for ''Green'' data'
81           READ (*,*) GP
82           IF (((RP*HP) .GT. NH) .OR. ((GP*HP) .GT. NH)) GOTO 150
83
84        ELSE
85           CALL GETFIL(IRUN,NT,NH,H,C,DXTIME,NFRAMS,
86     +          TITLE,LRG,.FALSE.,beam)
87        ENDIF
88
89        CALL FORBCK(C,NH,NT,F(1,1,K),B(1,1,K),RP,GP,HP)
90        X(K)=H
91  10  CONTINUE
92
93      WRITE(*,*)
94      CALL LOGQYN(' >  Another set of runs ?','N',LMORE)
95      IF (LMORE) GOTO 1
96      NRUNS=K
97
98      WRITE(*,*)
99   3  WRITE(*,120)
100 120  FORMAT(' >  First and last time channel ?  : ',$)
101      READ(*,*,ERR=3) I1,I2
102      IF (I1.LT.1 .OR. I2.GT.NT .OR. I2.LT.I1) GOTO 3
103      DO 20 K=1,NRUNS
104        CALL INTSYM(F(1,1,K),B(1,1,K),I1,I2,YI(K),EI(K))
105        CALL DIFSYM(F(1,1,K),B(1,1,K),I1,I2,YD(K),ED(K))
106  20  CONTINUE
107
108      WRITE(TITLE(21:39),130)I1,I2
109 130  FORMAT('  Bins: ',I4' - ',I4)
110      WRITE(TITLE(40:60),140)ISTART,IEND
111 140  FORMAT('  Runs: ',I5' - ',I5)
112
113      WRITE(*,*)
114      WRITE(*,*)' Integral Asymmetry'
115      CALL PLOT(X,YI,EI,NRUNS,1,SIZE,IWIDTH,TITLE)
116      WRITE(*,*)
117      WRITE(*,*)' Differential Asymmetry'
118      CALL PLOT(X,YD,ED,NRUNS,2,SIZE,IWIDTH,TITLE)
119      WRITE(*,*)
120      CALL LOGQYN(' >  Try another time window ?','Y',AGAIN)
121      WRITE(*,*)
122      IF (AGAIN) GOTO 3
123      END
124C
125      SUBROUTINE FORBCK(C,NH,NT,F,B,RP,GP,HP)
126C     ------------------------------
127C
128      INTEGER RP,GP,FB
129      REAL HP,C(NH,NT),F(2,*),B(2,*)
130C
131      FB=HP/2
132      DO 30 J=1,NT
133        F(1,J)=0.
134        B(1,J)=0.
135        F(2,J)=0.
136        B(2,J)=0.
137        do 29 K=1,FB
138          F(1,J)=F(1,J)+C(K+((RP-1)*HP),J)
139          F(2,J)=F(2,J)+C(K+((GP-1)*HP),J)
140          B(1,J)=B(1,J)+C(K+(((RP-1)*HP)+FB),J)
141          B(2,J)=B(2,J)+C(K+(((GP-1)*HP)+FB),J)
142  29    continue
143  30  CONTINUE
144      END
145C
146      SUBROUTINE INTSYM(F,B,I1,I2,YI,EI)
147C     ----------------------------------
148C
149      REAL F(2,*),B(2,*)
150C
151      SUMFR=0.0
152      SUMBR=0.0
153      SUMFG=0.0
154      SUMBG=0.0
155      DO 10 I=I1,I2
156        SUMFR=SUMFR+F(1,I)
157        SUMBR=SUMBR+B(1,I)
158        SUMFG=SUMFG+F(2,I)
159        SUMBG=SUMBG+B(2,I)
160  10  CONTINUE
161      YIF=(SUMFG-SUMFR)/(SUMFG+SUMFR)
162      YIB=(SUMBG-SUMBR)/(SUMBG+SUMBR)
163      YI=YIB-YIF
164      VARIF=(1.0+YIF*YIF)/(SUMFG+SUMFR)
165      VARIB=(1.0+YIB*YIB)/(SUMBG+SUMBR)
166      EI=SQRT(VARIF+VARIB)             ! VARIB prev YARIB
167      END
168C
169      SUBROUTINE DIFSYM(F,B,I1,I2,YD,ED)
170C     ----------------------------------
171C
172      REAL F(2,*),B(2,*)
173C
174      SUM=0.0
175      VARD=0.0
176      DO 10 I=I1,I2
177        FNORM=1.0/(F(2,I)+F(1,I))
178        ZF=(F(2,I)-F(1,I))*FNORM
179        BNORM=1.0/(B(2,I)+B(1,I))
180        ZB=(B(2,I)-B(1,I))*BNORM
181        SUM=SUM+ZB-ZF
182        VARD=VARD+(1.0+ZF*ZF)*FNORM+(1.0+ZB*ZB)*BNORM
183  10  CONTINUE
184      YD=SUM/FLOAT(I2-I1+1)
185      ED=SQRT(VARD)/FLOAT(I2-I1+1)
186      END
187C
188C***<plottting routines>************************************************
189C
190      SUBROUTINE ASKSIZ(SIZE,IWIDTH)
191C     ------------------------------
192C
193   1  WRITE(*,100)
194 100  FORMAT(' PLOT>  Character Size ? (def=1.5)  : ',$)
195      READ(*,200,ERR=1) NNN,SIZE
196 200  FORMAT(Q,F)
197      IF (NNN.EQ.0) SIZE=1.5
198      IF (SIZE.LT.0.001) GOTO 1
199   2  WRITE(*,110)
200 110  FORMAT(' PLOT>  Line-widths ?    (def=1; try 5 for PS) : ',$)
201      READ(*,210,ERR=2) NNN,IWIDTH
202 210  FORMAT(Q,I)
203      IF (NNN.EQ.0) IWIDTH=1
204      IF (IWIDTH.LT.1 .OR. IWIDTH.GT.21) GOTO 2
205      END
206C
207      SUBROUTINE PLOT(X,Y,E,N,ISYM,SIZE,IWIDTH,TITLE)
208C     -----------------------------------------------
209C
210      REAL      X(*),Y(*),E(*),YLOW(200),YHIGH(200)
211      LOGICAL   LPRINT,LWRITE,LERROR
212      CHARACTER FILNAM*60,TITLE*60
213C
214      CALL MAXMIN(X,Y,N,XMIN,XMAX,YMIN,YMAX)
215      WRITE(*,*)
216      CALL LOGQYN(' Plot>  Display the error-bars ?','N',LERROR)
217      WRITE(*,*)
218      CALL PGBEGIN(0,'?',1,1)
219      CALL PGSCH(SIZE)
220      CALL PGSLW(IWIDTH)
221      CALL PGENV(XMIN,XMAX,YMIN,YMAX,0,0)
222      CALL PGPOINT(N,X,Y,4)
223      IF (LERROR) THEN
224        DO 10 I=1,N
225          YLOW(I)=Y(I)-E(I)
226          YHIGH(I)=Y(I)+E(I)
227  10    CONTINUE
228        CALL PGERRY(N,X,YLOW,YHIGH,0.0)
229      ENDIF
230      IF (ISYM.EQ.1) THEN
231        CALL PGLABEL('Field  (G)','Integral Asymmetry',TITLE)
232      ELSEIF (ISYM.EQ.2) THEN
233        CALL PGLABEL('Field  (G)','Average Asymmetry',TITLE)
234      ENDIF
235      CALL PGEND
236
237      WRITE(*,*)
238      CALL LOGQYN(' Output>  Write out an ASCII file ?','N',LWRITE)
239      IF (LWRITE) THEN
240        WRITE(*,*)
241   1    WRITE(*,100)
242 100    FORMAT(' Output>  Filename ?  : ',$)
243        READ(*,200,ERR=1) FILNAM
244 200    FORMAT(A)
245        OPEN(UNIT=3,FILE=FILNAM,STATUS='NEW',FORM='FORMATTED',ERR=1)
246        DO 20 I=1,N
247          WRITE(3,110) X(I),Y(I),E(I)
248 110      FORMAT(X,F7.1,2F9.5)
249  20    CONTINUE
250        CLOSE(UNIT=3)
251      ENDIF
252
253      WRITE (*,*)
254      CALL LOGQYN(' Print>  Create a PGPLOT.PS file ?','N',LPRINT)
255      IF (LPRINT) THEN
256         CALL PGBEGIN(0,'/PS',1,1)
257         CALL PGSCH(SIZE)
258         CALL PGSLW(IWIDTH)
259         CALL PGENV(XMIN,XMAX,YMIN,YMAX,0,0)
260         CALL PGPOINT(N,X,Y,4)
261         IF (LERROR) THEN
262            DO 30 I=1,N
263               YLOW(I)=Y(I)-E(I)
264               YHIGH(I)=Y(I)+E(I)
265  30        CONTINUE
266            CALL PGERRY(N,X,YLOW,YHIGH,0.0)
267         ENDIF
268         IF (ISYM.EQ.1) THEN
269            CALL PGLABEL('Field  (G)','Integral Asymmetry',TITLE)
270         ELSEIF (ISYM.EQ.2) THEN
271            CALL PGLABEL('Field  (G)','Average Asymmetry',TITLE)
272         ENDIF
273         CALL PGEND
274      ENDIF
275
276      END
277C
278      SUBROUTINE MAXMIN(X,Y,N,XMIN,XMAX,YMIN,YMAX)
279C     ------------------------------------------
280C
281      REAL X(*),Y(*)
282C
283      XMIN=+1.0E20
284      XMAX=-1.0E20
285      YMIN=+1.0E20
286      YMAX=-1.0E20
287      DO 10 I=1,N
288        IF (X(I).LT.XMIN) XMIN=X(I)
289        IF (X(I).GT.XMAX) XMAX=X(I)
290        IF (Y(I).LT.YMIN) YMIN=Y(I)
291        IF (Y(I).GT.YMAX) YMAX=Y(I)
292  10  CONTINUE
293      XDIF=(XMAX-XMIN)/20.0+1.0E-10
294      YDIF=(YMAX-YMIN)/20.0+1.0E-10
295      XMIN=XMIN-XDIF
296      XMAX=XMAX+XDIF
297      YMIN=YMIN-YDIF
298      YMAX=YMAX+YDIF
299      WRITE(*,*)
300   1  WRITE(*,100) YMIN,YMAX
301 100  FORMAT(' Plot>  Ymin & Ymax ?  (def=',F7.3',',F7.3')  : ',$)
302      READ(*,200,ERR=1) NNN,YMIN1,YMAX1
303 200  FORMAT(Q,2F)
304      IF (NNN.NE.0) THEN
305        YMIN=YMIN1
306        YMAX=YMAX1
307      ENDIF
308      END
309C
310C***<logical question>**************************************************
311C
312      SUBROUTINE LOGQYN(S,D,L)
313C     ------------------------
314C
315      LOGICAL      L,L2
316      CHARACTER*1  D,D2,A
317      CHARACTER*45 STRING
318      CHARACTER    S(*)
319C
320      IF (D.EQ.'Y') THEN
321        L=.TRUE.
322        D2='N'
323        L2=.FALSE.
324      ELSEIF (D.EQ.'N') THEN
325        L=.FALSE.
326        D2='Y'
327        L2=.TRUE.
328      ELSE
329        WRITE(*,*)' Default should be Y or N !'
330        RETURN
331      ENDIF
332      CALL STPACK(STRING,S,45)
333   1  WRITE(*,100)STRING,D
334      READ(*,110)N,A
335      IF (N.EQ.0) RETURN
336      IF (A.EQ.'y' .OR. A.EQ.'T' .OR. A.EQ.'t') A='Y'
337      IF (A.EQ.'n' .OR. A.EQ.'F' .OR .A.EQ.'f') A='N'
338      IF (A.EQ.D) THEN
339        RETURN
340      ELSEIF (A.EQ.D2) THEN
341        L=L2
342        RETURN
343      ENDIF
344      GOTO 1
345 100  FORMAT(A45' Y/N (def=',A1')  : ',$)
346 110  FORMAT(Q,:,A1)
347      END
348
349C
350      SUBROUTINE STPACK(STRING,S,N)
351C     -----------------------------
352C
353      CHARACTER STRING(*),S(*)
354C
355      DO 10 I=1,N
356        STRING(I)=S(I)
357        IF (S(I).EQ.'?') GOTO 20
358  10  CONTINUE
359  20  DO 30 J=I+1,N
360  30    STRING(J)=' '
361      END
362C
363C***<hacked from Steve Johnston's program mcs_read.for>*****************
364C
365      SUBROUTINE GETFIL(IRUN,NT,NH,H,C,DXTIME,NFRAMS,
366     +                  TITLE,LRG,PRNT,beam)
367C     ---------------------------------------------------------
368C
369        include 'muon_def.inc'
370        integer*4 file
371        character*80 mcs_file
372        integer*4 rdata(128,2048)
373c
374        REAL      C(*),H
375        CHARACTER TITLE*60
376        LOGICAL   PRNT
377        integer lrg, beam
378C
379        file=irun
380
381c
382c  ***  manufacture file name and search VAX & HUB ***
383c
384        if(beam .eq. 1) then
385        write(mcs_file,1070)file
386        else if(beam .eq. 2) then
387        write(mcs_file,1071)file
388        else if(beam .eq. 3) then
389        write(mcs_file,1072)file
390        endif
391        goto 205
392201     if(beam .eq. 1) then
393        write(mcs_file,1060)file
394        else if(beam .eq. 2) then
395        write(mcs_file,1061)file
396        else if(beam .eq. 3) then
397        write(mcs_file,1062)file
398        endif
399        type *,'File not present in FEM, searching restore directory'
400        goto 205
401202     if(beam .eq. 1) then
402        write(mcs_file,1080)file
403        else if(beam .eq. 2) then
404        write(mcs_file,1081)file
405        else if(beam .eq. 3) then
406        write(mcs_file,1082)file
407        endif
408        type *,'File not present in restore directory, searching...'
409        goto 205
410203     type *,
411     + 'File not present. May need to be restored from the archive.'
412        type *,' '
413        STOP
414
415205     status=NXMread(mcs_file)
416
417        do ih=1,NXM_histogram_number
418           do ic=1,NXM_histogram_length
419              rdata(ih,ic) = NXM_histogram_counts(ih,ic)
420           enddo
421        enddo
422       
423        NH=NXM_histogram_number
424        NT=NXM_histogram_length
425        DXTIME=NXM_histogram_Resolution/1000.0
426        NFRAMS=NXM_beam_frames
427        LRG= NXM_run_switching_states
428        K=0
429        DO J=1,NT
430          DO I=1,NH
431            K=K+1
432            C(K)=FLOAT(RDATA(I,J))
433          ENDDO
434        ENDDO
435
436        if (PRNT) THEN
437c
438c  ***  print title and header of file  ***
439c
440           type *,' '
441           type *,' '
442           type *,'Sample                : ', NXM_sample_name
443           type *,'Temperature           : ', NXM_sample_temperature
444        endif
445
446        type *,'Field                 : ', NXM_sample_magnetic_field
447        H = NXM_sample_magnetic_field
448
449        if (PRNT) THEN
450           type *,'Orientation           : ', NXM_detector_orientation
451           type *,'Start time            : ', NXM_run_start_time
452           type *,'Stop time             : ', NXM_run_stop_time
453           type *,' '
454           type 1020, NXM_histogram_resolution
455           type 1030, NXM_beam_frames
456           type 1040, NXM_histogram_number
457           type 1050, NXM_histogram_length
458        endif
459
4601020    format(' Resolution            : ',1f10.2,' pS')
4611030    format(' Total no. of frames   : ',1i8)
4621040    format(' Number of histograms  : ',1i4)
4631050    format(' Number of time bins   : ',1i5)
464
465c       for EMU:
4661060    format('scratch$disk:[emumgr.restore]emu',i8.8,'.nxs') ! 2nd
4671070    format('emu$disk0:[data.emu]emu',i8.8,'.nxs') ! 1st
4681080    format('inst$disk:[emumgr.data]emu',i8.8,'.nxs') ! 3rd
469c       for MUSR:
4701061    format('scratch$disk:[mumgr.restore]musr',i8.8,'.nxs') ! 2nd
4711071    format('musr$disk0:[data.musr]musr',i8.8,'.nxs') ! 1st
4721081    format('inst$disk:[mumgr.data]musr',i8.8,'.nxs') ! 3rd
473c       for DEVA:
4741062    format('scratch$disk:[mutmgr.restore]mut',i8.8,'.nxs') ! 2nd
4751072    format('mut$disk0:[data.mut]mut',i8.8,'.nxs') ! 1st
4761082    format('inst$disk:[mutmgr.data]mut',i8.8,'.nxs') ! 3rd
477        end
478C
479CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC