Ticket #562: asym.for

File asym.for, 12.9 KB (added by Nick Draper, 11 years ago)
Line 
1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2C
3c (subroutines LOGQYN and STPACK copied from ASYM_RF.FOR 4/1/96.  PJCK)
4
5c 11/1/1999 J.S. Lord   Added DEVA file access
6
7      PROGRAM ASYM
8C     ------------
9C
10      REAL    X(200),YI(200),YD(200),C(65536),F(2048,200),B(2048,200)
11      REAL    DELTA(32),EI(200),ED(200)
12      LOGICAL AGAIN
13      integer inst
14      DATA    NCAL /32/
15C
16      CALL ASKSIZ(SIZE,IWIDTH)
17      CALL CALIB(DELTA,NCAL)
18      WRITE(*,*)
19c      CALL LOGQYN(' Instrument>  Are data from MuSR (else EMU) ?','Y',
20c     *             LMUSR)
21
2257      write(*,*)'Select instrument:'
23        write(*,*)'1. MuSR'
24        write(*,*)'2. EMU'
25        write(*,*)'3. DEVA'
26        read(*,*)inst
27        if(inst.lt.1 .or. inst.gt.3) goto 57
28
29      WRITE(*,*)
30   1  WRITE(*,100)
31 100  FORMAT(' >  First and last run number ?  : ',$)
32      READ(*,*,ERR=1) ISTART,IEND
33      IF (ISTART.LT.1 .OR. IEND.LT.ISTART) GOTO 1
34      NRUNS=IEND-ISTART+1
35      IF (NRUNS.GT.200) GOTO 1
36   2  WRITE(*,110)
37 110  FORMAT(' >  Balance parameter alpha ?    : ',$)
38      READ(*,*,ERR=2) ALPHA
39      IF (ALPHA.LT.0.01 .OR. ALPHA.GT.100.0) GOTO 2
40      K=1
41      DO 10 IRUN=ISTART,IEND
42        CALL GETFIL(IRUN,NT,NH,H,C,DXTIME,NFRAMS,inst)
43        IF (NCAL.LE.32 .AND. NCAL.NE.NH) STOP ' # of Hist mismatch!'
44        CALL FORBCK(C,NH,NT,ALPHA,F(1,K),B(1,K),DELTA,DXTIME,NFRAMS)
45        X(K)=H
46        K=K+1
47  10  CONTINUE
48      WRITE(*,*)
49   3  WRITE(*,120)
50 120  FORMAT(' >  First and last time channel ?  : ',$)
51      READ(*,*,ERR=3) I1,I2
52      IF (I1.LT.1 .OR. I2.GT.NT .OR. I2.LT.I1) GOTO 3
53      DO 20 K=1,NRUNS
54        CALL INTSYM(F(1,K),B(1,K),I1,I2,YI(K),EI(K))
55        CALL DIFSYM(F(1,K),B(1,K),I1,I2,YD(K),ED(K))
56  20  CONTINUE
57      WRITE(*,*)
58      WRITE(*,*)' Integral Asymmetry'
59      CALL PLOT(X,YI,EI,NRUNS,1,SIZE,IWIDTH)
60      WRITE(*,*)
61      WRITE(*,*)' Differential Asymmetry'
62      CALL PLOT(X,YD,ED,NRUNS,2,SIZE,IWIDTH)
63      WRITE(*,*)
64      CALL LOGQYN(' >  Try another time window ?','Y',AGAIN)
65      WRITE(*,*)
66      IF (AGAIN) GOTO 3
67      END
68C
69      SUBROUTINE CALIB(DELTA,NCAL)
70C     ----------------------------
71C
72      REAL         DELTA(*)
73      CHARACTER*60 FILNAM
74C
75      WRITE(*,*)
76   1  WRITE(*,100)
77 100  FORMAT(' INPUT>  Calibration Filename ?   : ',$)
78      READ(*,200,ERR=1) NNN,FILNAM
79 200  FORMAT(Q,A)
80      IF (NNN.EQ.0) THEN
81        WRITE(*,*)
82   2    WRITE(*,110)
83 110    FORMAT(' Calibration>  Dead time ? (nS)  : ',$)
84        READ(*,*,ERR=2) DTIME
85        DO 10 I=1,NCAL
86  10      DELTA(I)=DTIME
87        NCAL=1000
88      ELSE
89      OPEN(UNIT=2,FILE=FILNAM,STATUS='OLD',FORM='FORMATTED',READONLY,
90     *     ERR=2)
91        READ(2,*)
92        READ(2,*)
93        READ(2,*)
94        READ(2,*)
95        DO 20 I=1,NCAL
96  20      READ(2,*,ERR=3,END=3) DELTA(I)
97   3    NCAL=I-1
98        CLOSE(UNIT=2)
99        WRITE(*,*)
100        WRITE(*,*)' No. of detectors in calibration file = ',NCAL
101      ENDIF
102      END
103C
104      SUBROUTINE FORBCK(C,NH,NT,ALPHA,F,B,DELTA,DX,NFRAMS)
105C     ----------------------------------------------------
106C
107      REAL C(NH,NT),F(*),B(*),DELTA(*),DELTA2(32)
108C
109      DO 10 I=1,NH
110  10    DELTA2(I)=DELTA(I)/(DX*FLOAT(NFRAMS))
111      NH2=NH/2
112      DO 30 J=1,NT
113        SUMF=0.0
114        SUMB=0.0
115        DO 20 I=1,NH2
116          IB=I+NH2
117          DI=C(I,J)
118          DIB=C(IB,J)
119          SUMF=SUMF+DI/(1.0-DI*DELTA2(I))
120          SUMB=SUMB+DIB/(1.0-DIB*DELTA2(I))
121  20    CONTINUE
122        F(J)=SUMF
123        B(J)=ALPHA*SUMB
124  30  CONTINUE
125      END
126C
127      SUBROUTINE INTSYM(F,B,I1,I2,YI,EI)
128C     ----------------------------------
129C
130      REAL F(*),B(*)
131C
132      SUMF=0.0
133      SUMB=0.0
134      DO 10 I=I1,I2
135        SUMF=SUMF+F(I)
136        SUMB=SUMB+B(I)
137  10  CONTINUE
138      YI=(SUMF-SUMB)/(SUMF+SUMB)
139      VARI=(1.0+YI*YI)/(SUMF+SUMB)
140      EI=SQRT(VARI)
141      END
142C
143      SUBROUTINE DIFSYM(F,B,I1,I2,YD,ED)
144C     ----------------------------------
145C
146      REAL F(*),B(*)
147C
148      SUM=0.0
149      VARD=0.0
150      DO 10 I=I1,I2
151        FBNORM=1.0/(F(I)+B(I))
152        Z=(F(I)-B(I))*FBNORM
153        SUM=SUM+Z
154        VARD=VARD+(1.0+Z*Z)*FBNORM
155  10  CONTINUE
156      YD=SUM/FLOAT(I2-I1+1)
157      ED=SQRT(VARD)/FLOAT(I2-I1+1)
158      END
159C
160C***<plottting routines>************************************************
161C
162      SUBROUTINE ASKSIZ(SIZE,IWIDTH)
163C     ------------------------------
164C
165   1  WRITE(*,100)
166 100  FORMAT(' PLOT>  Character Size ? (def=1.5)  : ',$)
167      READ(*,200,ERR=1) NNN,SIZE
168 200  FORMAT(Q,F)
169      IF (NNN.EQ.0) SIZE=1.5
170      IF (SIZE.LT.0.001) GOTO 1
171   2  WRITE(*,110)
172 110  FORMAT(' PLOT>  Line-widths ?    (def=1; try 5 for PS) : ',$)
173      READ(*,210,ERR=2) NNN,IWIDTH
174 210  FORMAT(Q,I)
175      IF (NNN.EQ.0) IWIDTH=1
176      IF (IWIDTH.LT.1 .OR. IWIDTH.GT.21) GOTO 2
177      END
178C
179      SUBROUTINE PLOT(X,Y,E,N,ISYM,SIZE,IWIDTH)
180C     -----------------------------------------
181C
182      REAL         X(*),Y(*),E(*),YLOW(200),YHIGH(200)
183      LOGICAL      LWRITE,LERROR
184      CHARACTER*60 FILNAM
185C
186      CALL MAXMIN(X,Y,N,XMIN,XMAX,YMIN,YMAX)
187      WRITE(*,*)
188      CALL LOGQYN(' Plot>  Display the error-bars ?','N',LERROR)
189      WRITE(*,*)
190      CALL PGBEGIN(0,'?',1,1)
191      CALL PGSCH(SIZE)
192      CALL PGSLW(IWIDTH)
193      CALL PGENV(XMIN,XMAX,YMIN,YMAX,0,0)
194      CALL PGPOINT(N,X,Y,4)
195      IF (LERROR) THEN
196        DO 10 I=1,N
197          YLOW(I)=Y(I)-E(I)
198          YHIGH(I)=Y(I)+E(I)
199  10    CONTINUE
200        CALL PGERRY(N,X,YLOW,YHIGH,0.0)
201      ENDIF
202      IF (ISYM.EQ.1) THEN
203        CALL PGLABEL('Field  (G)','Asymmetry','Integral Asymmetry')
204      ELSEIF (ISYM.EQ.2) THEN
205        CALL PGLABEL('Field  (G)','Asymmetry','Differential Asymmetry')
206      ENDIF
207      CALL PGEND
208      WRITE(*,*)
209      CALL LOGQYN(' Output>  Write out an ASCII file ?','N',LWRITE)
210      IF (LWRITE) THEN
211        WRITE(*,*)
212   1    WRITE(*,100)
213 100    FORMAT(' Output>  Filename ?  : ',$)
214        READ(*,200,ERR=1) FILNAM
215 200    FORMAT(A)
216        OPEN(UNIT=3,FILE=FILNAM,STATUS='NEW',FORM='FORMATTED',ERR=1)
217        DO 20 I=1,N
218          WRITE(3,110) X(I),Y(I),E(I)
219 110      FORMAT(X,F7.1,2F9.5)
220  20    CONTINUE
221        CLOSE(UNIT=3)
222      ENDIF
223      END
224C
225      SUBROUTINE MAXMIN(X,Y,N,XMIN,XMAX,YMIN,YMAX)
226C     ------------------------------------------
227C
228      REAL X(*),Y(*)
229C
230      XMIN=+1.0E20
231      XMAX=-1.0E20
232      YMIN=+1.0E20
233      YMAX=-1.0E20
234      DO 10 I=1,N
235        IF (X(I).LT.XMIN) XMIN=X(I)
236        IF (X(I).GT.XMAX) XMAX=X(I)
237        IF (Y(I).LT.YMIN) YMIN=Y(I)
238        IF (Y(I).GT.YMAX) YMAX=Y(I)
239  10  CONTINUE
240      XDIF=(XMAX-XMIN)/20.0+1.0E-10
241      YDIF=(YMAX-YMIN)/20.0+1.0E-10
242      XMIN=XMIN-XDIF
243      XMAX=XMAX+XDIF
244      YMIN=YMIN-YDIF
245      YMAX=YMAX+YDIF
246      WRITE(*,*)
247   1  WRITE(*,100) YMIN,YMAX
248 100  FORMAT(' Plot>  Ymin & Ymax ?  (def=',F7.3',',F7.3')  : ',$)
249      READ(*,200,ERR=1) NNN,YMIN1,YMAX1
250 200  FORMAT(Q,2F)
251      IF (NNN.NE.0) THEN
252        YMIN=YMIN1
253        YMAX=YMAX1
254      ENDIF
255      END
256CC
257C***<logical question>**************************************************
258C
259      SUBROUTINE LOGQYN(S,D,L)
260C     ------------------------
261C
262      LOGICAL      L,L2
263      CHARACTER*1  D,D2,A
264      CHARACTER*45 STRING
265      CHARACTER    S(*)
266C
267      IF (D.EQ.'Y') THEN
268        L=.TRUE.
269        D2='N'
270        L2=.FALSE.
271      ELSEIF (D.EQ.'N') THEN
272        L=.FALSE.
273        D2='Y'
274        L2=.TRUE.
275      ELSE
276        WRITE(*,*)' Default should be Y or N !'
277        RETURN
278      ENDIF
279      CALL STPACK(STRING,S,45)
280   1  WRITE(*,100)STRING,D
281      READ(*,110)N,A
282      IF (N.EQ.0) RETURN
283      IF (A.EQ.'y' .OR. A.EQ.'T' .OR. A.EQ.'t') A='Y'
284      IF (A.EQ.'n' .OR. A.EQ.'F' .OR .A.EQ.'f') A='N'
285      IF (A.EQ.D) THEN
286        RETURN
287      ELSEIF (A.EQ.D2) THEN
288        L=L2
289        RETURN
290      ENDIF
291      GOTO 1
292 100  FORMAT(A45' Y/N (def=',A1')  : ',$)
293 110  FORMAT(Q,:,A1)
294      END
295
296C
297      SUBROUTINE STPACK(STRING,S,N)
298C     -----------------------------
299C
300      CHARACTER STRING(*),S(*)
301C
302      DO 10 I=1,N
303        STRING(I)=S(I)
304        IF (S(I).EQ.'?') GOTO 20
305  10  CONTINUE
306  20  DO 30 J=I+1,N
307  30    STRING(J)=' '
308      END
309C
310C***<hacked from Steve Johnston's program mcs_read.for>*****************
311C
312        SUBROUTINE GETFIL(IRUN,NT,NH,H,C,DXTIME,NFRAMS,inst)
313C       -----------------------------------------------------
314C
315        REAL    C(*)
316        integer inst
317C
318        character*80 mcs_file
319        integer*4 file
320        integer*4 ncan
321        integer*2 ib(512)
322        integer*2 istfla,rescod,nrun,numhis,histlen,numrec,reclen,i
323        integer*2 t0bin(32),tgood_begin(32),tgood_end(32),j
324        integer*2 ngroups,branch(4),crate(4),first(4),ndev(4)
325        integer*2 begin,end
326        integer*4 res_pscnds,ntotal(32)
327        character rtitle*40, devtyp(4)*4
328        character*9 start_date,stop_date
329        character*8 start_time,stop_time
330        logical*1 rgmode
331        integer*2 nreads, framread
332        character*32 grouping
333        integer*4 rdata
334        dimension rdata(64,2048)
335        integer a,rtot,asy
336        real redgrn,ntot(32),alpha,totred,totgrn,temp
337        character*60 comment
338c
339        equivalence
340        1 (istfla,ib(1)), (rescod,ib(2)), (res_pscnds,ib(7)),
341        1 (nrun,ib(4)), (start_date,ib(110)), (stop_date,ib(115)),
342        1 (start_time,ib(120)), (stop_time,ib(124)), (rtitle,ib(70)),
343        1 (numhis,ib(16)), (histlen,ib(15)), (numrec,ib(65)),
344        1 (reclen,ib(66)), (t0bin(1),ib(256-3*32)),
345        1 (tgood_begin(1),ib(256-2*32)), (tgood_end(1),ib(256-32)),
346        1 (ntotal(1),ib(256)), (branch(1),ib(512-5*4)),
347        1 (crate(1),ib(512-4*4)), (first(1),ib(512-3*4)),
348        1 (ndev(1),ib(512-2*4)), (devtyp(1),ib(512-4)),
349        1 (ngroups,ib(512)),(comment,ib(431)),
350        1 (rgmode,ib(67)),(nreads,ib(68)),(framread,ib(69))
351C
352        file=0
353        type *,' '
354c
355c  ***  manufacture file name and search VAX & HUB ***
356c
357        FILE=irun
358        type *,' '
359        redgrn=1
360        IF (inst.eq.1) THEN
361          write(mcs_file,1070)file
362        ELSE IF (inst.eq.2) THEN
363          write(mcs_file,1071)file
364        ELSE
365          write(mcs_file,1072)file
366        ENDIF
367        open(1,file=mcs_file,status='OLD',form='UNFORMATTED',
368     1  err=201,readonly)
369        goto 205
370201     IF (inst.eq.1) THEN
371          write(mcs_file,1060)file
372        ELSE IF (inst.eq.2) THEN
373          write(mcs_file,1061)file
374        ELSE
375          write(mcs_file,1062)file
376        ENDIF
377        type *,'File not present in uVAX,searching the HUB'
378        open(1,file=mcs_file,status='OLD',form='UNFORMATTED',
379     1  err=202,readonly)
380        goto 205
381202     IF (inst.eq.1) THEN
382          write(mcs_file,1080)file
383        ELSE IF (inst.eq.2) THEN
384          write(mcs_file,1081)file
385        ELSE
386          write(mcs_file,1082)file
387        ENDIF
388        type *,'File not present in uVAX,searching the HUB'
389        open(1,file=mcs_file,status='OLD',form='UNFORMATTED',
390     1  err=203,readonly)
391        goto 205
392203     type *,'File not present in HUB,re-enter valid file number'
393        type *,' '
394        RETURN
395c
396c  ***  enter data from MCS file into decimal form (from binary)  ***
397c
398205     read(1) ib
399        if (rgmode.eq.0) goto 99
400        redgrn=2
40199      do 100 ih=1,numhis*redgrn
402100        read(1) (rdata(ih,i),i=1,histlen)
403        close(1)
404        rtitle(30:30)=rtitle(29:29)
405c
406c  ***  print title and header of file  ***
407c
408        type *,' '
409102     rtitle(11:11)=' '
410        rtitle(20:20)='K'
411        rtitle(29:29)='G'
412        type *,' '
413        type *,'Sample                : ',rtitle(1:11)
414        type *,'Temperature           : ',rtitle(12:20)
415        type *,'Field                 : ',rtitle(21:29)
416        READ(RTITLE(21:28),*) H
417        if (rtitle(30:30).ne.'l') goto 801
418        type *,'Orientation           : LONGITUDINAL'
419        goto 803
420801     if (rtitle(30:30).ne.'t') goto 802
421        type *,'Orientation           : TRANSVERSE'
422        goto 803
423802     type *,'Orientation           : ',rtitle(30:30)
424        numgroup=5
425803     type *,'Collimation           : ',rtitle(31:40)
426        type *,' '
427        type *,'Start time            : ',start_time
428        type *,'Stop time             : ',stop_time
429        type *,'Start date            : ',start_date
430        type *,'Stop date             : ',stop_date
431        type *,' '
432        type 1020,res_pscnds
433        ncan=int(nreads)*int(framread)
434        type 1030,ncan
435        type 1040,numhis
436        type 1050,histlen
437        NH=numhis
438        NT=histlen
439        DXTIME=FLOAT(res_pscnds)/1000.0
440        NFRAMS=ncan
441        K=0
442        DO J=1,NT
443          DO I=1,NH
444            K=K+1
445            C(K)=FLOAT(RDATA(I,J))
446          ENDDO
447        ENDDO
4481020    format(' Resolution            : ',1i6,' pS')
4491030    format(' Total no. of frames   : ',1i8)
4501040    format(' Number of histograms  : ',1i3)
4511050    format(' Number of time bins   : ',1i5)
4521060    format('MUSR$DISK:[MUMGR.DATA]MUS',i5.5,'.RAW')
4531070    format('MUSR$DISK0:[DATA.MUSR]R',i5.5,'.RAL')
4541080    format('SCRATCH$DISK:[MUMGR.RESTORE]USR',i5.5,'.RAW')
4551061    format('MUSR$DISK:[EMUMGR.DATA]EMU',i5.5,'.RAW')
4561071    format('EMU$DISK0:[DATA.EMU]R',i5.5,'.RAL')
4571081    format('SCRATCH$DISK:[EMUMGR.RESTORE]EMU',i5.5,'.RAW')
4581062    format('MUT$DISK:[MUTMGR.DATA]MUT',i5.5,'.RAW')
4591072    format('MUT$DISK0:[DATA.MUT]R',i5.5,'.RAL')
4601082    format('SCRATCH$DISK:[MUTMGR.RESTORE]MUT',i5.5,'.RAW')
461        end
462C
463CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC