Ticket #4513: TFITD.FOR;74

File TFITD.FOR;74, 43.8 KB (added by Peter Parker, 9 years ago)

Tfit Routine

Line 
1*
2* Fit tof spectrum to sum of N Voigt or Gaussian functions.
3* Fitting parameters are heights and sigmas in y. Position fixed by mass.
4* Individual widths or areas can also be fixed
5      PROGRAM TFIT_FSE
6      PARAMETER (NP=20,ndet=200,npts=1024)
7      REAL CM(NP,NP),CMI(NP,NP),ERR(NP),DE1G(3),DE1L(3)
8      REAL WORK(23000),F(3000),P(30),PS(30),PST(30)
9      REAL TH(200),RL0(200),RL1(200),DT0A(200),
10     $M,K0,K1,L0,L1,SIGY(20)
11      REAL COMP(20),DCOMP(20),WIDW(20),DWIDW(20),AV(20),DAV(20)
12     $,RV(20),DRV(20)
13      REAL xva(ndet,npts),yva(ndet,npts),eva(ndet,npts)
14      real xin(npts),yin(npts),ein(npts)
15      real xout(npts),yout(npts),eout(npts),yout2(np,npts)
16      REAL wt(np,ndet),wte(np,ndet),sig(np,ndet)
17     *,sige(np,ndet),Chi(ndet)
18   
19      COMMON QA(3000,3),E0A(3000,3),YA(20,3000,3),TMAX
20     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
21      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
22      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
23     $ DTE(20,3),TR0(20),NPEAKS
24      CHARACTER MSDATA*70, FNAMEIN*80
25      CHARACTER DET1*1,DET2*2,DET3*3,FNAME*80,fin*70,FDAT*80
26
27C  READ IN DATA.
28      write(6,*) ' Name of file containing time of flight data?'
29      read(5,'(a)') FIN
30      write(6,5) fin
31    5 format(' input file is ',80a)
32      FDAT=FIN//'D.DAT'
33    6 write(6,*) ' tmin,tmax? (tmax-tmin must be 2**n)'
34      read(5,*) tmin,tmax
35      npt=tmax-tmin
36      if(npt.ne.64.and.npt.ne.128.and.npt.ne.256.and.npt.ne.512)then
37       write(6,*) ' tmin-tmax must be 2**n'
38       go to 6
39      end if
40      ndmin=1000
41      ndmax=0
42      open(unit=3,file=fin,status='old')
43  10  read(3,*) idet,lpt
44*      write(6,*) ' idet=',idet,' lpt=',lpt
45      if(idet.gt.1000) go to 20
46*      WRITE(6,*) ' IDET=',IDET
47      if(idet.lt.ndmin)ndmin=idet
48      if(idet.gt.ndmax)ndmax=idet
49      iv=0
50      do i=1,lpt
51       read(3,*) xv,yv,ev
52       if(yv.eq.0.0) yv=1e-6
53       if(ev.eq.0.0) ev=1e-6
54*       write(6,*) xv,yv,ev
55       if(xv.ge.tmin.and.xv.lt.tmax)then
56        iv=iv+1
57        xva(idet,iv)=xv
58        yva(idet,iv)=yv
59        eva(idet,iv)=ev
60*       write(6,*) iv,xv,yv,ev
61       end if
62      end do
63      lpt=iv
64      lptlast=lpt
65      go to 10
66   20 irun1=idet
67      write(6,*) ' irun1=',irun1
68      lpt=lptlast
69      close(3)
70      write(6,*) ' number of points=',lpt
71      write(6,*) ' ndmin=',ndmin,' ndmax=',ndmax
72
73      write(6,*) ' Name of file containing fit parameters?'
74      read(5,'(a)') fnamein
75C  TAKE INPUT FROM PREVIOUSLY PREPARED DATA FILE
76      OPEN(UNIT=3,FILE=fnamein,STATUS='OLD')
77      E1(1)=4897.4 ! Calibrated value of E1
78
79* DEFINE ENERGY RESOLUTION FUNCTION AND IPNO
80* OPTION 1 IS FORWARD SCATTERING
81* OPTION 2 IS BACK SCATTERING SINGLE DIFFERENCE, DETECTORS 3-90.
82* OPTION 3 IS BACK SCATTERING DOUBLE DIFFERENCE* DETECTORS 3-90.
83* OPTION 4 IS BACK SCATTERING SINGLE DIFFERENCE, DETECTORS 91-135.
84* OPTION 5 IS BACK SCATTERING DOUBLE DIFFERENCE* DETECTORS 91-135.
85*
86      READ(3,*) IOPT
87      IF(IOPT.EQ.1) THEN
88      DE1G(1)=74.0
89      DE1L(1)=24.0
90*      write(6,*) ' de1g=',de1g(1),' de1l=',de1l(1)
91      IF(irun1.lt.15210)IPNO=0004
92      IF(irun1.gt.15210)IPNO=0005
93      ELSE IF(IOPT.EQ.5) THEN
94      DE1G(1)=54.0
95      DE1L(1)=137.0
96      IF(irun1.lt.15210)IPNO=0004
97      IF(irun1.gt.15210)IPNO=0005
98      ELSE IF(IOPT.EQ.3) THEN
99      DE1G(1)=54.0
100      DE1L(1)=137.0
101      IF(irun1.lt.15210)IPNO=0004
102      IF(irun1.gt.15210)IPNO=0005
103      ELSE IF(IOPT.EQ.4) THEN
104      DE1G(1)=54.0
105      DE1L(1)=137.0
106      IF(irun1.lt.15210)IPNO=0004
107      IF(irun1.gt.15210)IPNO=0005
108      ELSE IF(IOPT.EQ.2) THEN
109      DE1G(1)=88.0
110      DE1L(1)=41.0
111      IF(irun1.lt.15210)IPNO=0004
112      IF(irun1.gt.15210)IPNO=0005
113      END IF
114      write(6,*) ' irun1=',irun1,' ipno=',ipno
115
116*      write(6,*) ' e1=',e1(1)
117      E1(2)=60300.0
118      E1(3)=78400.0
119     
120*      READ(3,*) DE1G(1),DE1L(1) ! Energy resolution.
121*      write(6,*) ' de1g=',de1g(1),' de1l=',de1l(1)
122      DE1L(2)=202.0
123      DE1L(3)=192.0
124      DE1G(2)=1.0
125      DE1G(3)=1.0
126*      READ(3,*) IPNO ! Instrument parameter file number.
127*      write(6,*) ' ipno=',ipno
128*      READ(3,*) MAXFUN
129*      write(6,*) ' maxfun=',maxfun
130      MAXFUN=1000
131      READ(3,*) NPEAKS ! Number of different masses in sample.
132      write(6,*) ' npeaks=',npeaks
133      DO I=1,NPEAKS
134      READ(3,*) AM(I),XS(I),WID(I)
135      write(6,*) i,' am=',am(i),' xs=',xs(i),' wid=',wid(i)
136      END DO
137
138*      READ(3,*) NBACK
139      NBACK=0
140*      write(6,*) ' nback=',nback
141      IF(NBACK.EQ.1) THEN
142       WRITE(6,*) ' BACKGROUND SUBTRACTED'
143      ELSE
144*       WRITE(6,*) ' NO BACKGROUND SUBTRACTED'
145      END IF
146
147      NBPTS=0
148*      READ(3,*) NBPTS
149      IF(NBPTS.EQ.1) THEN
150        WRITE(6,*) ' ERRORS ON BAD PTS SET TO 1E6'
151      ELSE
152*        WRITE(6,*) ' BAD PTS NOT EXCLUDED'
153      END IF
154
155*      READ(3,*) NMSSUB
156*       write(6,*) ' nmssub=',nmssub
157      NMSSUB=0
158      IF(NMSSUB.EQ.1) THEN
159       WRITE(6,*) ' MULTIPLE SCATTERING SUBTRACTED'
160       READ(3,'(70A)') MSDATA
161       write(6,'(70a)') msdata
162       READ(3,*) NFACT
163       IF(NFACT.EQ.1) THEN
164        WRITE(6,*) ' MS MUTIPLIED BY FITTING FACTOR'
165       END IF
166       CALL MSNORM(XI,YI,LPT,MSDATA,ND,RMS)
167      ELSE
168*       WRITE(6,*) ' NO MULIPLE SCATTERING SUBTRACTED'
169      END IF
170      READ(3,*) NFSE
171      IF(NFSE.EQ.1) THEN
172       WRITE(6,*) ' FSE SUBTRACTED'
173      ELSE
174       WRITE(6,*) ' FSE NOT SUBTRACTED'
175      END IF
176      CLOSE(3)
177
178
179
180
181   25 write(6,*) ' First and last spectrum number?'
182      read(5,*) nfirst,nlast
183      if(nfirst.lt.ndmin.or.nlast.gt.ndmax) then
184       write(6,*) ' first spectrum must be >',ndmin
185       write(6,*) ' last spectrum must be <',ndmax
186       go to 25
187      end if
188     
189      OPEN(UNIT=2,FILE=FDAT,STATUS='NEW')
190      do nvv=nfirst,nlast   ! Start of loop over spectra.
191*      open(unit=4,file='temp.dat',status='new')       
192       n=nvv
193      do it=1,lpt
194       xin(it)=xva(n,it)
195       yin(it)=yva(n,it)
196       ein(it)=eva(n,it)
197*       write(4,*) xin(it),yin(it),ein(it)
198      end do
199*      close(4)
200       
201      TCHW=XIN(10)-XIN(9)
202      WRITE(6,*) ' TIME CHANNEL WIDTH=',TCHW
203      WRITE(6,*) ' NUMBER OF INPUT DATA POINTS=',LPT
204      YMAX=0.0
205      sum=0.0
206      DO I=1,LPT
207       IF(YIN(I).GT.YMAX) YMAX=YIN(I)
208      END DO
209      DO I=1,LPT
210       IF(YIN(I).EQ.0.0) THEN
211        YIN(I)=YMAX*1E-6
212        EIN(I)=1E6
213       END IF
214      END DO
215     
216      DO I=1,LPT
217      XI(I)=XIN(I)
218      YI(I)=YIN(I)
219      EI(I)=EIN(I)
220      END DO
221
222      ND=N
223      NPL=1
224      PI=ACOS(-1.0)
225
226C 2.  READ IN VALUES OF THETA AND L1 FROM FILE
227      CALL PAR_READ(TH,DT0A,RL0,RL1,IPNO,NDET,NS)
228     
229
230C 3. DEFINE FIXED RESOLUTION COMPONENTS
231      CALL RES_READ(DL0,DL1,DT0,DW)
232      DL0=DL0/100.0
233      DL1=DL1/100.0
234      DW=DW/100.0
235      DO I=1,NS
236      TH(I)=TH(I)*PI/180.0
237      END DO
238
239      ISPEC=N
240      L0=RL0(N)
241      L1=RL1(N)
242      THETA=TH(N)
243      DTHETA=0.5*DW/RL1(N)
244      THETAW=TH(N)*180.0/PI
245      DTHW=DTHETA*180.0/PI
246      T0=DT0A(N)
247      if(iopt.eq.2.and.ispec.ge.3.and.ispec.le.46) then
248       de1g(1)=93.0
249       de1l(1)=33.0
250      else if(iopt.eq.2.and.ispec.ge.47.and.ispec.le.90) then
251       de1g(1)=88.0
252       de1l(1)=41.0
253      else if(iopt.eq.2.and.ispec.ge.91.and.ispec.le.134) then
254       de1g(1)=85.0
255       de1l(1)=47.5
256      end if
257
258      WRITE(6,*) ' DETECTOR NUMBER=',N
259      WRITE(6,*) ' L0=',L0,' +-',DL0,' METRES'
260      WRITE(6,*) ' L1=',L1,' +-',DL1,' METRES'
261      WRITE(6,*) ' THETA=',THETAW,' +-',DTHW,' DEGREES'     
262      WRITE(6,*) ' E1=',E1
263      WRITE(6,*) ' DE1L=',DE1L(1),' DE1G=',DE1G(1)
264      WRITE(6,*) ' T0=',T0,' +-',DT0, 'uSEC'
265      WRITE(6,*) ' NUMBER OF PEAKS=',NPEAKS
266*      WRITE(6,*) ' MAXFUN=',MAXFUN
267*      write(6,*) l0,l1,theta,e1,t0,am(i)
268
269
270      NXSP=0 ! Number of peak amplitude fitting parameters
271      DO I=1,NPEAKS
272      IF(XS(I).EQ.0.0) NXSP=NXSP+1
273      IF(AM(I).EQ.0.0) GO TO 130
274      END DO
275      IF(NXSP.EQ.NPEAKS) XS(1)=1.0 ! amp parameters not independent.
276
277*  Calculate tof at at peak centre.
278*      write(6,*) l0,l1,theta,e1,t0,am(i)
279      DO I=1,NPEAKS
280      TR0(I)=TREC(L0,L1,THETA,E1,T0,AM(I))
281      END DO
282     
283* 5. Calculate Q and E0 and y for each tof and each energy.
284
285      DO IE=1,1
286       V1=SQRT(E1(IE)/5.2276E-6)
287       T1=L1/V1
288      DO I=1,LPT
289*       T=XIN(I)+TCHW ! Fudge
290       T=XIN(I)
291       V0=L0/(T*1E-6-T0*1E-6-T1)
292       E0=5.2276E-6*V0**2
293       K0=SQRT(E0/2.0717)
294       K1=SQRT(E1(IE)/2.0717)
295       QA(I,IE)=SQRT(K0**2+K1**2-2.0*K0*K1*COS(THETA))
296       E0A(I,IE)=E0
297      END DO
298      END DO
299
300* Calculate y for each t and for each mass.
301      DO IE=1,1
302      DO IM=1,NPEAKS
303      DO I=1,LPT
304      W=E0A(I,IE)-E1(IE)
305      RM=AM(IM)/1.00867 ! Mass in multiples of neutron mass.
306      WR=2.0717*QA(I,IE)**2/RM
307      YA(IM,I,IE)=0.2413*(RM/QA(I,IE))*(W-WR)
308      END DO
309      END DO
310      END DO
311
312      QAV=QA(LPT/2,1)
313
314* 6. Calculate resolution function for each mass in time of flight.
315      DO IE=1,1
316      DO IM=1,NPEAKS
317       DTG(IM,IE)= TRESN(L0,DL0,L1,DL1,THETA,DTHETA
318     $ ,E1(IE),DE1G(IE),T0,DT0,AM(IM))
319       DTE(IM,IE)= TRESN(L0,0.0,L1,0.0,THETA,0.0,
320     $ E1(IE),DE1L(IE),T0,0.0,AM(IM))
321      END DO   
322      END DO
323      DT=XIN(2)-XIN(1)
324      DO IE=1,1
325      DO IM=1,NPEAKS
326       DO IT=1,LPT/2
327        TV=(IT-1)*DT
328        TRES(IM,IT,IE)=DT*VOIGT(TV,DTG(IM,IE),DTE(IM,IE),0)
329       END DO
330       DO IT=LPT/2+1,LPT
331        TV=(IT-1-LPT)*DT
332        TRES(IM,IT,IE)=DT*VOIGT(TV,DTG(IM,IE),DTE(IM,IE),0)
333       END DO
334      END DO
335      END DO
336
337      WRITE(6,*) ' '
338      WRITE(6,*) '   MASS     ','  POSITION ','    INTENSITY ',
339     $'   WIDTH    ','    DTG     ','    DTE     '
340      DO I=1,NPEAKS
341      WRITE(6,50) AM(I),TR0(I),XS(I),WID(I),DTG(I,1),DTE(I,1)
342      if(tr0(i).lt.tmin.or.tr0(i).gt.tmax) then
343       write(6,*) ' Peak for mass', Am(i),' lies outside fitting range'
344       stop
345      end if
346   50 FORMAT(' ',1P6E12.4)
347      END DO
348
349C Calculate estimates for amplitude.
350      YMAX=0.0
351      DO I=1,LPT
352      IF(YMAX.LT.YI(I)) THEN
353      YMAX=YI(I)
354      IMAX=I
355      END IF
356      END DO
357
358* Estimate of standard deviation in tof.
359      DO I=1,NPEAKS
360      WD=800/11.604 ! Estimated Debye energy in meV.
361      SIGY(I)=SQRT(0.75*0.1196*AM(I)*WD)
362      END DO
363
364      NPAR=0 ! Number of fitting parameters other than scale and background.
365      DO I=1,NPEAKS
366*       write(6,*) i, ' xs=',xs(i),' wid=',wid(i)
367       IF(XS(I).EQ.0.0) THEN
368        NPAR=NPAR+1
369        P(NPAR)=10.0*ymax*AM(I)/SIGY(I)
370        WRITE(6,*) ' PARAMETER',NPAR,' IS AMP OF PEAK',I
371*        WRITE(6,*) ' Start value=',p(npar)
372       END IF
373       IF(WID(I).EQ.0.0) THEN
374        NPAR=NPAR+1
375        P(NPAR)=SIGY(I)
376        WRITE(6,*) ' PARAMETER', NPAR, ' IS WIDTH OF PEAK',I
377       END IF
378      END DO
379      N=NPAR+1
380      P(N)=YMAX ! Scale Factor.
381      NSC=N ! Number of parameter containing scale factor.
382      WRITE(6,*) ' PARAMETER',N,' IS SCALE FACTOR'
383      CALL CTS(YOUT,LPT,NPEAKS,P)
384
385
386      NV=0
387      nr=0
388      DO I=1,NPEAKS
389      SUMD=0.0
390      SUMF=0.0
391*       WRITE(6,*) I, ' M=',AM(I),' TR=',tr0(i)
392       DO IT=1,LPT
393*        write(6,*) it,xin(it),yin(it),yout(it)
394         trp=tr0(I)+2.0*TCHW
395         trm=tr0(i)-2.0*TCHW
396*        write(6,*) it,' t=',xin(it),' trm=',trm,' trp=',trp 
397        IF(XIN(IT).GT.TRM.AND.xin(it).LT.TRP) THEN
398         SUMD=SUMD+YIN(IT)
399         SUMF=SUMF+YOUT(IT)
400*         write(6,*) ' t=',xin(it),' data=',yin(it),' estimate=',yout(it) 
401*        write(6,*) it,' t=',xin(it),' trm=',trm,' trp=',trp 
402*        write(6,*) it,'sumd=',sumd,' sumf=',sumf
403        END IF
404       END DO
405       if(sumf.eq.0.0) then
406        write(6,*) ' Bad Spectrum'
407        do it=1,lpt
408         yout(it)=0.0
409         eout(it)=0.0
410        end do
411        go to 100
412       end if
413
414       RATIO=SUMD/SUMF
415       IF(XS(I).EQ.0.0) THEN
416        NV=NV+1
417        P(NV)=P(NV)*RATIO
418        WRITE(6,*) ' PARAMETER',nv,' =',p(nv),' IS AMP OF PEAK',I
419*        write(6,*) ' SUMD=',SUMD,' SUMF=',SUMF,' ratio=',ratio
420       ELSE
421        nr=nr+1
422        if(nr.lt.2)P(N)=P(N)*RATIO
423       END IF
424       IF(WID(I).EQ.0.0) THEN
425        nv=nv+1
426        WRITE(6,*) ' PARAMETER',nv,' =',p(nv),' IS width OF PEAK',I
427       END IF
428      END DO
429
430      IF(NBACK.EQ.1) THEN
431       N=N+3
432       CALL BSTART(XIN,YIN,LPT,TMAX,A,B,C)
433       P(N-2)=A
434       P(N-1)=B
435       P(N)=C
436      END IF
437      IF(NMSSUB.EQ.1.AND.NFACT.EQ.1) THEN
438      N=N+1
439      P(N)=1.0
440      WRITE(6,*) ' PARAMETER',N,' IS MS FACTOR'
441      END IF
442
443      WRITE(6,*) ' Number of fitting parameters=',N
444
445
446C SCALE PARAMETERS
447      DO 70 I=1,N
448      V(I)=1.0/P(I)
449      PS(I)=P(I)*V(I)
450      PST(I)=P(I)! Store start values of fit parameters.
451  70  CONTINUE
452
453C  CALCULATE DMAX USED IN VA05A.
454      DMAX=3.0
455C  H IS DISTANCE BETWEEN X PTS USED TO CALCULATE PARTIAL DERIVATIVES.
456      H=DMAX*1E-2
457C
458      ACC=0.001
459      NDMAX=0
460   80 IPRINT=0
461
462      CALL VA05A(LPT,N,F,PS,H,DMAX,ACC,MAXFUN,IPRINT,WORK)
463
464      IF(IPRINT.EQ.1.AND.NDMAX.NE.1) THEN
465      TYPE *,' SUM OF SQUARES FAILED TO DECREASE'
466      TYPE *, ' DMAX REDUCED TO 0.1'
467      write(6,*) ' ndmax=',ndmax
468      NDMAX=1
469      DMAX=0.1
470      H=DMAX*1E-2
471      write(6,*) ' ndmax=',ndmax
472      GO TO 80
473      END IF
474
475      IF(IPRINT.EQ.2) TYPE *,' MAXIMUM NUMBER OF CALLS TO VAO5A MADE '
476
477      LINC=0
478      CS=CHISQC(PS,N,LPT,LINC) ! Calculate chi-square
479      CS=CS/(LINC-N) ! Reduced chi-square.
480      WRITE(6,*) ' FIT COMPLETED CHI=',CS
481
482      WRITE(6,*) 'ERRORS NOW BEING CALCULATED.'
483
484      if(maxfun.ne.0.and.cs.lt.100.0) then
485* Calculate errors on fitted parameters.
486      CALL ERRORS(N,LPT,PS,CM,CMI,ERR)
487
488* Check calculated errors.
489      CALL ERRCHECK(N,LPT,PS,ERR)
490      else
491      do i=1,n
492       err(i)=1e6
493      end do
494      end if
495
496      DO 90 I=1,N
497      P(I)=PS(I)/V(I)
498      ERR(I)=ERR(I)/V(I)
499   90 CONTINUE
500
501      DO I=1,N
502      WRITE(6,*) I, ' PST=',PST(I),' P=',P(I),' +-',err(i)
503      END DO
504
505      NP1=0 ! Number of peak amplitude parameters
506      NPARV=0 ! Number of fitting parameters
507      SUMF=0.0 ! Sum of fixed amplitudes.
508      DO I=1,NPEAKS
509       IF(XS(I).EQ.0.0) THEN
510        NP1=NP1+1
511        NPARV=NPARV+1
512        AV(NP1)=ABS(P(NPARV))
513        DAV(NP1)=ABS(ERR(NPARV))
514       ELSE
515        SUMF=SUMF+ABS(XS(I))
516       END IF
517       IF(WID(I).EQ.0.0) NPARV=NPARV+1
518      END DO
519      NPF=NP1+1 ! Number of amplitude fitting parameters.
520      AV(NPF)=SUMF*ABS(P(NSC))
521      DAV(NPF)=SUMF*ABS(ERR(NSC))
522
523
524      SUM=0.0
525      DO I=1,NPF
526       SUM=SUM+ABS(AV(I))
527      END DO
528
529      DO M=1,NPF
530       SUM1=0.0
531       SUM2=0.0
532       DO N=1,NPF
533        IF(N.NE.M) THEN
534         SUM1=SUM1+AV(N)
535         SUM2=SUM2+DAV(N)**2
536        END IF
537       END DO
538       DRSQ=(DAV(M)**2*SUM1**2+AV(M)**2*SUM2)/sum**4
539       RV(M)=AV(M)/SUM
540       DRV(M)=SQRT(DRSQ)
541      END DO
542
543      NP1=0
544      NPARV=0
545      DO I=1,NPEAKS
546       IF(XS(I).EQ.0.0) THEN
547        NP1=NP1+1
548        NPARV=NPARV+1
549        COMP(I)=RV(NP1)
550        DCOMP(I)=DRV(NP1)
551       ELSE
552        COMP(I)=ABS(XS(I)*RV(NPF))/SUMF
553        DCOMP(I)=ABS(XS(I)*DRV(NPF))/SUMF
554       END IF
555       
556       IF(WID(I).EQ.0.0) THEN
557        NPARV=NPARV+1
558        WIDW(I)=P(NPARV)
559        DWIDW(I)=ERR(NPARV)
560       ELSE
561        WIDW(I)=WID(I)
562        DWIDW(I)=1E-6
563       END IF
564       WIDW(I)=ABS(WIDW(I)) ! Eliminate negative values
565       DWIDW(I)=ABS(DWIDW(I)) ! Eliminate negative values
566       
567      END DO
568
569      write(6,*) ' scale factor=',p(npar+1)
570
571  100 Continue
572      WRITE(6,*) '       MASS    ','  RELATIVE WT ',
573     $'   ERROR    ',' SD OF J(Y) ','   ERROR', '       CHISQ'
574      WRITE(6,*) ' spectrum number=',ND
575      DO I=1,NPEAKS
576      WRITE(6,95) I,AM(I),COMP(I),DCOMP(I),WIDW(I),DWIDW(I),CS
577   95 FORMAT(' ',I3,1P6E12.4)
578      wt(i,nd)=comp(i)
579      wte(i,nd)=dcomp(i)
580      sig(i,nd)=widw(i)
581      sige(i,nd)=dwidw(i)
582      end do
583      Chi(nd)=CS
584   
585     
586      PMIN=0.00001
587      IF(NBPTS.EQ.1) THEN
588       CALL CTS(YOUT,LPT,NPEAKS,P)
589       DO I=1,LPT
590        WDIFF=ABS((YIN(I)-YOUT(I))/EIN(I))
591        PROB=2.0*(1-ERF(WDIFF))
592        IF(PROB.LT.PMIN) THEN
593         EIN(I)=1E6
594        WRITE(6,*) ' T=',XIN(I),WDIFF,PROB
595        END IF
596       END DO
597      END IF
598
599      CALL CTSF(YOUT,LPT,NPEAKS,P)
600      CALL CTSF2(YOUT2,LPT,NPEAKS,P)
601   
602     
603      write(2,*) NVV,lpt,npeaks
604      DO I=1,LPT
605       write(2,*) xin(i),yin(i),EIN(I),YOUT(I),(yout2(ip,i),ip=1,npeaks) 
606      END DO
607   97 format(' ',1e12.4,6e11.3)
608      IF(NMSSUB.EQ.1.AND.NFACT.EQ.1) THEN
609       IF(ND.LT.10) THEN
610        WRITE(DET1(1:1),'(I1.1)')ND
611        FNAME='TEMP'//DET1
612       ELSE IF(ND.GE.10.AND.ND.LT.100) THEN
613        WRITE(DET2(1:2),'(I2.2)')ND
614        FNAME='TEMP'//DET2
615       ELSE
616        WRITE(DET3(1:3),'(I3.3)')ND
617        FNAME='TEMP'//DET3
618       END IF
619     
620*       WRITE(6,2) FNAME
621*    2  FORMAT(' ',' MS FACTOR IN FILE;',80A)
622*       OPEN(UNIT=2,FILE=FNAME,STATUS='NEW')
623*       WRITE(2,*) P(N)
624*       CLOSE(2)
625      END IF
626  130 CONTINUE 
627
628      END DO ! END OF LOOP OVER SPECTRA
629      CLOSE(2)
630     
631* Create output file containing fitted parameters.
632      ndv=nlast-nfirst+1
633     
634      do i=nfirst,nlast
635*      write(6,*) i, wt(1,i),wte(1,i)
636      end do
637
638      CALL FOUT(NPEAKS,Nfirst,ndv,AM,WT,WTE,SIG,SIGE,CHI,FIN,IPNO)
639
640
641      STOP
642      END
643
644      SUBROUTINE CALFUN(M,N,F,PS)
645      REAL F(3000),PS(30),P(30)
646      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
647     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
648      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
649      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
650     $ DTE(20,3),TR0(20),NPEAKS
651
652      DO 10 I=1,N
653      P(I)=PS(I)/V(I)
654      if(p(i).lt.0.0)p(i)=-p(i)
655   10 CONTINUE     
656     
657      CHI=0.0
658      CALL CTS(F,M,NPEAKS,P)
659      DO 20 I=1,M
660      X=XI(I)
661      if(ei(i).ne.0.0) F(I)=(F(I)-YI(I))/EI(I)
662      CHI=CHI+F(I)**2
663   20 CONTINUE
664      CHI=CHI/M
665   30 FORMAT(7E10.3)
666
667      RETURN
668      END
669
670C Function to calculate tof at recoil peak in micro seconds
671C L0 = incident flight path in metres.
672C L1= final flight path in metres.
673C TH = scattering angle in radians.
674C E1 = analyser energy in meV.
675C DT0 = time delay in microseconds.
676C M=atomic mass in amu.
677      FUNCTION TREC(L0,L1,TH,E1,DT0,M)
678      REAL L0,L1,M
679      M=M/1.00867 ! Convert to multiple of neutron mass.
680      ARG=M**2-SIN(TH)**2
681      IF(ARG.LT.0.0) ARG=0.0
682      FACT=(COS(TH)+SQRT(ARG))/(M+1)
683      M=M*1.00867
684      V1=SQRT(E1/5.2276E-6)
685      V0=V1/FACT
686      TREC=L0/V0+L1/V1+DT0*1E-6
687      TREC=TREC*1E6
688      RETURN
689      END
690 
691C Function to calculate y.
692C L0 = incident flight path in metres.
693C L1= final flight path in metres.
694C TH = scattering angle in radians.
695C T = time of flight in microsec.
696C E1 = analyser energy in meV.
697C DT0 = time delay in microseconds.
698C M=atomic mass in amu.
699       FUNCTION Y(L0,L1,TH,E1,T,DT0,M)
700       REAL L0,L1,M,K0,K1
701       M=M/1.00867 ! Convert to multiple of neutron mass.
702       V1=SQRT(E1/5.2276E-6)
703       T1=L1/V1
704       T0=T*1E-6-DT0*1E-6-T1
705       V0=L0/T0
706       E0=5.2276E-6*V0**2
707       W=E0-E1
708       K0=SQRT(E0/2.0717)
709       K1=SQRT(E1/2.0717)
710       Q=SQRT(K0**2+K1**2-2.0*K0*K1*COS(TH))
711       WR=2.0717*Q**2/M
712       Y=0.2413*M*(W-WR)/Q
713       M=M*1.00867
714       RETURN
715       END
716
717C
718C    Voigt function centred at X=X0
719C    sigma is gaussian standard deviation
720c    DYE is Lorentzian DYE.
721C    Peak area normalised to 1.
722C    WRITTEN BY WIFD
723C    Modified by JM
724C
725      FUNCTION VOIGT(X,SIGMA,DYE,X0)
726
727      DOUBLE PRECISION WR,WI,XX,YY
728C
729     
730      GAMMA=DYE*2.0
731      XS=X
732      X=XS-X0
733      OVRTPI=0.564189584
734      OVRT2=0.707106781
735      BTEM=OVRT2/SIGMA
736      ATEM=OVRTPI*BTEM
737      XTEM=X*BTEM
738      YTEM=0.5*GAMMA*BTEM
739      XX= DBLE(XTEM)
740      YY= DBLE(YTEM)
741      CALL WERF(WR,WI,XX,YY)
742      SWR=SNGL(WR)
743      SWI=SNGL(WI)
744      CTEM=ATEM*BTEM
745      VOIGT=ATEM*SWR
746      X=XS
747      DWRDX=-2.*(XTEM*SWR-YTEM*SWI)
748      DWRDY= 2.*(YTEM*SWR+XTEM*SWI-OVRTPI)
749      DERX=CTEM*DWRDX
750      DERS=-ATEM*(SWR+DWRDX*XTEM+DWRDY*YTEM)/SIGMA
751      DERG=0.5*CTEM*DWRDY
752C
753      RETURN
754      END
755c
756c
757        subroutine WERF(rs1,rs2,xx,yy)
758c       W.I.F.David 25-May-84
759        implicit real*8         (a-h,o-z)
760        real*8                  lambda
761        logical                 b
762        x=dabs(xx)
763        y=dabs(yy)
764        if (y .lt. 4.29 .and. x .lt. 5.33) go to 1
765        h= 0.
766        nc= 0
767        nu= 8
768        lambda= 0.
769        b= .true.
770        go to 2
771 1      s=(1.0-y/4.29)*dsqrt(1.0-x**2/28.41)
772        h=1.6*s
773        h2=2.0*h
774        nc=6+idint(23.0*s)
775        nu=9+idint(21.0*s)
776        lambda=h2**nc
777        b= .false.
778        if (lambda .eq. 0.) b= .true.
779 2      r1=0.
780        r2=0.
781        s1=0.
782        s2=0.
783        n=nu+1
784 3      n=n-1
785        fn=n+1
786        t1=y+h+fn*r1
787        t2=x-fn*r2
788        c=0.5/(t1**2+t2**2)
789        r1=c*t1
790        r2=c*t2
791        if (h .le. 0.0 .or. n .gt. nc) go to 4
792        t1= lambda+s1
793        s1=r1*t1-r2*s2
794        s2=r2*t1+r1*s2
795        lambda=lambda/h2
796 4      if (n .gt. 0) go to 3
797        if (b) go to 6
798        rs1=s1
799        rs2=s2
800        go to 7
801 6      rs1=r1
802        rs2=r2
803 7      rs1= 1.12837916709551*rs1
804        if (y .eq. 0.0) rs1= dexp(-x**2)
805        rs2= 1.12837916709551*rs2
806        if (xx .lt. 0) rs2= -rs2
807        return
808        end
809
810
811      FUNCTION CHISQ(PS,N,M)
812     
813      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
814     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
815      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
816      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
817     $ DTE(20,3),TR0(20),NPEAKS
818      REAL PS(30),P(30),F(3000)
819
820     
821      DO I=1,N
822      P(I)=PS(I)/V(I)
823      END DO
824
825      CHISQ=0.0
826      CALL CTS(F,M,NPEAKS,P)
827      DO 20 I=1,M
828      X=XI(I)
829      IF(EI(I).ne.0.0) F(I)=(F(I)-YI(I))/EI(I)
830      IF(EI(I).LT.1E5) THEN ! Eliminate excluded points.
831       CHISQ=CHISQ+F(I)**2
832      END IF
833   20 CONTINUE
834
835      END
836
837      FUNCTION CHISQC(PS,N,M,LINC)
838     
839      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
840     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
841      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
842      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
843     $ DTE(20,3),TR0(20),NPEAKS
844      REAL PS(30),P(30),F(3000)
845
846      LINC=0 
847      DO I=1,N
848      P(I)=PS(I)/V(I)
849      END DO
850
851      CHISQC=0.0
852      CALL CTS(F,M,NPEAKS,P)
853      DO 20 I=1,M
854      X=XI(I)
855      IF(Ei(I).ne.0.0)F(I)=(F(I)-YI(I))/EI(I)
856      CHISQC=CHISQC+F(I)**2
857      IF(EI(I).LT.1E5)LINC=LINC+1
858   20 CONTINUE
859
860      END
861
862      SUBROUTINE ERRCHECK(N,M,PS,ERR)
863      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
864     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
865      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
866      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
867     $ DTE(20,3),TR0(20),NPEAKS
868      REAL PS(N),ERR(N)
869     
870      DO I=1,N
871      CS=CHISQ(PS,N,M)
872      PS(I)=PS(I)+ERR(I)
873      CSP=CHISQ(PS,N,M)
874      PS(I)=PS(I)-ERR(I)
875      DCS=CSP-CS
876      IF(DCS.LE.0.0) THEN
877      ERR(I)=PS(I)
878      WRITE(6,*) ' FOR PARAMETER',I,
879     $' CHISQ DECREASES AWAY FROM MINIMUM'
880      ELSE IF(DCS.LT.1.0) THEN
881      ERR(I)=ERR(I)/DCS
882      WRITE(6,*) ' PROBLEM WITH ERROR ON PARAMETER',I
883      WRITE(6,*) ' CHANGE IN CHISQ AT P(I)+DP(I)=',DCS
884      END IF
885      END DO
886
887      END
888
889      SUBROUTINE ERRORS(N,M,PS,CM,CMI,ERR)
890
891      REAL PS(N),ERR(N),CM(N,N),PD(30),CMI(N,N),DP(30)
892      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
893     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
894      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
895      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
896     $ DTE(20,3),TR0(20),NPEAKS
897
898      DO I=1,N
899      PD(I)=PS(I)
900      END DO
901
902* Find increments used in calculation of curvature matrix.
903      CALL INCREMENT(PS,N,M,DP)
904
905*  Calculate curvature matrix.
906      CALL CURVATURE(PS,DP,N,M,CM)
907
908*  Invert curvature matrix
909      CALL MINV(N,CM,CMI,SM)
910     
911*  CALCULATE ERRORS
912      DO I=1,N
913      IF(SM.EQ.0.0) ERR(I)=SQRT(ABS(CMI(I,I)))
914      IF(SM.EQ.1) ERR(I)=PS(I) ! Singular error matrix.
915      END DO
916
917      END
918
919*  Subroutine to invert matrix. Uses NAG subroutine F01AAF.
920*  A is N x N matrix. AI is inverse of A.
921      SUBROUTINE MINV(N,A,AI,SM)
922      REAL A(N,N),AI(N,N),W(50)
923
924* Check for singularity.
925      SM=0.0
926      DO I=1,N
927      ZERO=0.0
928      DO J=1,N
929      IF(A(I,J).NE.0.0) ZERO=1
930      END DO
931      IF(ZERO.EQ.0.0) THEN
932      WRITE(6,*) ' ALL ELEMENTS IN ROW',I,
933     $'OF ERROR MATRIX ARE ZERO'
934      SM=1.0
935      RETURN
936      END IF
937      END DO
938
939      IFAIL=-1
940      CALL F01AAE(A,N,N,AI,N,W,IFAIL)
941      IF(IFAIL.EQ.1) THEN
942      WRITE(6,*) ' ERROR MATRIX IS SINGULAR'
943      DO I=1,N
944      DO J=1,N
945      AI(I,J)=0.0
946      END DO
947      END DO
948      RETURN
949      END IF
950
951      END
952
953* Calculate increments for calculation of curvature
954* take increment which increases chisq by BETWEEN 1 AND 2.
955      SUBROUTINE INCREMENT(PS,N,M,DP)
956      REAL PS(N),DP(N),PD(30)
957      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
958     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
959      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
960      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
961     $ DTE(20,3),TR0(20),NPEAKS
962
963      DO I=1,N
964      PD(I)=PS(I)
965      END DO
966
967      CS=CHISQ(PS,N,M)
968
969      DO I=1,N
970
971*  Find value of DCS >2.
972      DP(I)=PS(I) ! Start value for increment
973      DPMIN =0
974      NINC=0
975   10 PD(I)=PS(I)+DP(I)
976
977      NINC=NINC+1 ! Check for no minimum.
978      IF(NINC.GT.10) THEN
979      DP(I)=PS(I)
980      GO TO 30
981      END IF
982     
983      CSP=CHISQ(PD,N,M)
984      DCS=CSP-CS ! Calculate change in chisq in step dp from min.
985      IF(DCS.LT.1.0) THEN
986       DPMIN=DP(I)
987       DP(I)=DP(I)*2.0
988       GO TO 10
989      ELSE IF(DCS.GT.3.0) THEN
990       DP(I)=DP(I)/2.0
991       go to 10
992      ELSE IF(DCS.GT.2.0) THEN
993       DPMAX=DP(I)
994      ELSE
995       GO TO 30
996      END IF
997
998* Now find value of dp(i) which changes chisq by between 1 and 2.
999   20 DP(I)=(DPMAX+DPMIN)/2.0
1000* Calculate change in chi-sq corresponding to dp(i).
1001      PD(I)=PS(I)+DP(I)
1002      CSP=CHISQ(PD,N,M)
1003      DCS=CSP-CS
1004      IF(DCS.GT.2.0) THEN
1005      DPMAX=DP(I)
1006      GO TO 20
1007      ELSE IF(DCS.LT.1.0) THEN
1008      DPMIN=DP(I)
1009      GO TO 20
1010      END IF
1011   30 CONTINUE
1012      PD(I)=PS(I) ! Reset pd to value at minimum.
1013
1014      END DO
1015      END
1016
1017
1018      SUBROUTINE CURVATURE(PS,DP,N,M,CM)
1019      REAL PS(N),DP(N),CM(N,N),PD(30)
1020      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
1021     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
1022      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
1023      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
1024     $ DTE(20,3),TR0(20),NPEAKS
1025
1026      DO I=1,N
1027      PD(I)=PS(I) ! Initialise pd.
1028      END DO
1029
1030      C=CHISQ(PS,N,M) ! chisq at minimum.
1031
1032      DO I=1,N
1033      DO J=1,N
1034      IF(J.LT.I) THEN
1035      CM(I,J)=CM(J,I)
1036      TIJ=0
1037      GO TO 5
1038      END IF
1039
1040      DPI=DP(I)
1041      DPJ=DP(J)
1042
1043      IF(I.NE.J) THEN
1044      PD(I)=PS(I)+DPI
1045      PD(J)=PS(J)+DPJ
1046      CIJ=CHISQ(PD,N,M)
1047      PD(I)=PS(I)
1048      PD(J)=PS(J)+DPJ
1049      CJ=CHISQ(PD,N,M)
1050      PD(I)=PS(I)+DPI
1051      PD(J)=PS(J)
1052      CI=CHISQ(PD,N,M)
1053      PD(I)=PS(I)
1054      PD(J)=PS(J)
1055      CM(I,J)=0.5*(CIJ-CI-CJ+C)/(DPI*DPJ)
1056      TIJ=(CIJ-C)/(DPI*DPJ)
1057
1058      ELSE IF(I.EQ.J) THEN
1059      PD(I)=PS(I)+2.0*DPI
1060      CII=CHISQ(PD,N,M)
1061      PD(I)=PS(I)+DPI
1062      CI=CHISQ(PD,N,M)
1063      CM(I,I)=0.5*(CII-2.0*CI+C)/(DPI**2)
1064      TIJ=(CII-C)/DPI**2
1065      END IF
1066    5 CONTINUE   
1067 
1068      END DO
1069      END DO
1070
1071      END
1072
1073* Reads in resolution parameters.
1074      SUBROUTINE RES_READ(DL0,DL1,DT0,DW)
1075      OPEN(UNIT=3,FILE='EVS$DISK0:[EVSMGR.userprogs]RESOLUTION.DAT',
1076     $STATUS='OLD',READONLY,SHARED)
1077      READ(3,*) DL0 ! Uncertainty in incident flight path (cm)
1078      READ(3,*) DL1 ! Uncertainty in scattered flight path(cm)
1079      READ(3,*) DT0 ! Tof uncertainty (usec)
1080      READ(3,*) DW  ! Detector width (cm)
1081      CLOSE(3)
1082      END
1083
1084* IPNO is IP run number,NDET is maximum no of detectors,NS is
1085* number of detectors in IPNRUN.dat.
1086      SUBROUTINE PAR_READ(TH,DT0,L0,L1,IPNO,NDET,NS)
1087      REAL TH(NDET),L0(NDET),L1(NDET),DT0(NDET)
1088      CHARACTER RUN*4,FIN*40
1089
1090* Define file name
1091      WRITE(RUN(1:4),'(I4.4)')IPNO
1092      FIN='EVS$disk0:[EVSMGR.CALIB.PAR]IP'//RUN
1093      WRITE(6,5) FIN
1094    5 FORMAT(' ',' IPFILE=',40A)
1095
1096      OPEN(UNIT=3,FILE=FIN,STATUS='OLD',READONLY,SHARED)
1097      NS=1
1098   10 READ(3,*,END=20) I,TH(NS),DT0(NS),L0(NS),L1(NS)
1099      NS=NS+1
1100      GO TO 10
1101   20 CLOSE(3)
1102      NS=NS-1
1103
1104      END
1105
1106
1107* Calculate time of flight spectrum.
1108* NM=number of masses, M=mass, XS=amplitude, SIGY=width in y.
1109* TRES=resolution fn in tof, L0,L1, lengths in metres, TH angle in rad.
1110* T0 time delay in usec, E1 analyser energy in meV.
1111      SUBROUTINE CTS(CT,NPTS,NPEAKS,P)
1112      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
1113     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
1114      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
1115      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
1116     $ DTE(20,3),TR0(20),NPEAKS1
1117      REAL CT(NPTS),M
1118      REAL RES(3000),CTM(3000),P(30),AMP(30),WIDTH(30)
1119
1120      NP=0
1121      DO I=1,NPEAKS
1122       IF(XS(I).EQ.0.0) THEN
1123        NP=NP+1
1124        AMP(I)=ABS(P(NP))
1125       ELSE
1126        AMP(I)=ABS(P(NPAR+1)*XS(I))
1127       END IF
1128       IF(WID(I).EQ.0.0)THEN
1129        NP=NP+1
1130        WIDTH(I)=P(NP)
1131       ELSE
1132        WIDTH(I)=WID(I)
1133       END IF
1134      END DO
1135
1136      DO IT=1,NPTS
1137       CT(IT)=0.0
1138      END DO
1139
1140      DO IE=1,1
1141       IF(IE.EQ.1) FACT=1.0
1142       IF(IE.EQ.2) FACT=0.516
1143       IF(IE.EQ.3) FACT=0.079
1144      DO IM=1,NPEAKS
1145       IF(AM(IM).LT.3) THEN ! Harmonic Oscillator
1146        DSQV=12.0*4.18036*WIDTH(IM)**4/AM(IM)
1147       ELSE ! Debye Approximation
1148        DSQV=12.8*4.18036*WIDTH(IM)**4/AM(IM)
1149       END IF
1150       if(nfse.eq.0) dsqv=0.0
1151       DO IT=1,NPTS
1152        YV=YA(IM,IT,IE)
1153        E0=E0A(IT,IE)
1154        Q=QA(IT,IE)
1155        IF(IE.EQ.1) THEN
1156         CPV=CP(YV,Q,AM(IM),WIDTH(IM),DSQV)
1157        ELSE
1158*         CPV=CP2(YV,yv2,Q,AM(IM),WIDTH(IM),DSQV)
1159         cpv=0.0 
1160        END IF
1161        CTM(IT)=AM(IM)*AMP(IM)*CPV*E0**0.1/Q
1162        RES(IT)=TRES(IM,IT,IE)
1163       END DO
1164
1165       CALL C06EKE(1,CTM,RES,NPTS,IFAIL) ! Convolve with resn fn.
1166       DO IT=1,NPTS
1167        CT(IT)=CT(IT)+FACT*CTM(IT)
1168       END DO
1169
1170      END DO
1171      END DO
1172
1173      IF(NBACK.EQ.1) THEN
1174       N=NP+4
1175       DO IT=1,NPTS
1176        X=XI(IT)
1177       CT(IT)=CT(IT)+P(N-2)*(TMAX-X)+P(N-1)*(TMAX-X)**2
1178     $ +P(N)*(TMAX-X)**3
1179       END DO                   
1180      END IF
1181
1182      IF(NMSSUB.EQ.1.AND.NFACT.EQ.1) THEN ! Muliple scattering correction.
1183       DO I=1,NPTS
1184        CT(I)=CT(I)+P(N+1)*RMS(I)
1185       END DO
1186      ELSE IF(NMSSUB.EQ.1.AND.NFACT.NE.1) THEN
1187      DO I=1,NPTS
1188       CT(I)=CT(I)+RMS(I)
1189      END DO
1190      END IF
1191
1192* Fudge to avoid unphysical widths in fit.       
1193       do im=1,npeaks
1194        IF(AM(IM).GT.7.0.AND.WIDTH(IM).LT.3.0) then
1195         do i=1,npts
1196          CT(I)=CT(I)*1.0e2
1197         END DO
1198        END IF
1199       END DO
1200
1201
1202      END
1203
1204
1205
1206
1207C Function to calculate resolution in time.
1208C L0 +- DL0 = incident flight path in metres.
1209C L1 +- DL1 = final flight path in metres.
1210C TH +- DTH = scattering angle in radians.
1211C E1 +- DE1 = analyser energy in meV.
1212C T0 +- DT0 = time delay in microseconds.
1213C M=atomic mass in amu.
1214* Calls TREC.
1215      FUNCTION TRESN(L0,DL0,L1,DL1,TH,DTH,E1,DE1,T0,DT0,M)
1216      REAL L0,L1,M
1217
1218      TR=TREC(L0,L1,TH,E1,T0,M)
1219*      write(6,*) ' tr=',tr,' t0=',t0
1220      DTL0=TREC(L0+DL0,L1,TH,E1,T0,M)-TR
1221      DTL1=TREC(L0,L1+DL1,TH,E1,T0,M)-TR
1222      DTTH=TREC(L0,L1,TH+DTH,E1,T0,M)-TR
1223      DTE1=TREC(L0,L1,TH,E1+DE1,T0,M)-TR
1224      DTT0=TREC(L0,L1,TH,E1,T0+DT0,M)-TR
1225      TRESN=SQRT(DTL0**2+DTL1**2+DTTH**2+DTE1**2+DTT0**2)
1226*      WRITE(6,*) ' L0=',L0,' +-',DL0,' L1=',L1,' +-',DL1
1227*      WRITE(6,*) ' TH=',TH,' +-',DTH,' E1=',E1,' =-',DE1
1228*      WRITE(6,*) ' T0=',T0,' +-',DT0, ' M=',M
1229*      WRITE(6,*) ' TRESN=',TRESN
1230
1231      END   
1232     
1233* Calculates J(y) including first Sears Correction terms
1234      FUNCTION CP2(Y1,y2,Q,M,SIG1,D2V)
1235      REAL Y,M,J0,J3
1236
1237      NP=101 ! Number of points in integration
1238      CP=0.0
1239      SIG=ABS(SIG1)
1240      DY=(y2-y1)/(np-1)
1241      do i=1,np
1242       y=y1+dy*(i-1)
1243       X=Y/SIG
1244       J0=EXP(-X**2/2)/(2.506628*SIG) ! 2.5066=SQRT(2PI)
1245       D3J=X*(3-X**2)*J0/SIG**3
1246       J3=-M*D2V*D3J/(150.49*Q)
1247       if(d2v.ne.0.0) then
1248        CP=CP+J0+J3
1249       else
1250        cp=cp+j0
1251       end if
1252      END DO
1253      CP2=CP*dy
1254
1255      END
1256
1257
1258
1259      SUBROUTINE BSTART(XIN,YIN,LPT,TMAX,A,B,C)
1260      REAL XIN(LPT),YIN(LPT),N1,N2
1261      T1=1E6
1262      T3=0.0
1263      Y1=0.0
1264      DO I=1,LPT
1265       IF(XIN(I).GT.497) THEN
1266        TMAX=XIN(I)
1267        YMAX=YIN(I)
1268       END IF
1269       IF(XIN(I).LT.T1) THEN
1270        DO J=I,10
1271         Y1=Y1+YIN(I)
1272        END DO
1273        T1=XIN(I)
1274        Y1=Y1/10.0
1275       END IF
1276       IF(XIN(I).LT.400.0.AND.XIN(I).GT.T3)THEN
1277        T3=XIN(I)
1278        Y3=YIN(I)
1279       END IF
1280      END DO
1281
1282      TMEAN=(T1+T3)/2
1283      T2=0.0
1284      DO I=1,LPT
1285       IF(XIN(I).GT.T2.AND.XIN(I).LE.TMEAN) THEN
1286        T2=XIN(I)
1287        Y2=YIN(I)
1288       END IF
1289      END DO
1290
1291*      WRITE(6,*) ' TMAX=',TMAX,' YMAX=',YMAX   
1292*      WRITE(6,*) ' T1=',T1,' Y1=',Y1 
1293*      WRITE(6,*) ' T2=',T2,' Y2=',Y2 
1294*      WRITE(6,*) ' T3=',T3,' Y3=',Y3 
1295
1296      X1=TMAX-T1
1297      X2=TMAX-T2
1298      X3=TMAX-T3
1299     
1300      D1=(Y1/X1**3-Y2/X2**3)/(1/X1-1/X2)
1301      D2=(Y1/X1**3-Y3/X3**3)/(1/X1-1/X3)
1302      N1=(1/X1**2-1/X2**2)/(1/X1-1/X2)
1303      N2=(1/X1**2-1/X3**2)/(1/X1-1/X3)
1304*      write(6,*) ' d1=',d1,' d2=',d2
1305*      write(6,*) ' n1=',n1,' n2=',n2
1306      A=(D1-D2)/(N1-N2)
1307*      write(6,*) ' a=',a
1308      B=((Y1/X1**3-Y2/X2**3)-A*(1/X1**2-1/X2**2))/(1/X1-1/X2)
1309*      write(6,*) ' b=',b
1310      C=Y1/X1**3-A/X1**2-B/X1
1311*      WRITE(6,*) ' A=',A,' B=',B,' C=',C
1312
1313      END
1314
1315
1316
1317      FUNCTION CP(Y,Q,M,SIG1,D2V)
1318      REAL Y,M,J0,J3
1319     
1320      SIG=ABS(SIG1)
1321      X=Y/SIG
1322      J0=EXP(-X**2/2)/(2.506628*SIG) ! 2.5066=SQRT(2PI)
1323      D3J=X*(3-X**2)*J0/SIG**3
1324      J3=-M*D2V*D3J/(150.49*Q)
1325      if(d2v.ne.0.0)then
1326      CP=J0+J3
1327      else
1328      CP=J0
1329      END IF
1330
1331      END
1332
1333
1334
1335      SUBROUTINE MSNORM(XI,YI,LPT,MSDATA,ND,RMS)
1336      CHARACTER MSDATA*70
1337      REAL XI(LPT),YI(LPT),RMS(LPT)
1338      real d(3000,20),de(3000,20),sum(20),X(3000),MSC(3000)
1339      REAL WORK1(3000),WORK2(3000),TOTSC(3000)
1340
1341      CHARACTER DET1*1,DET2*2,DET3*3,FNAME*80
1342
1343      DO I=1,3000
1344       MSC(I)=0.0
1345       TOTSC(I)=0.0
1346      END DO
1347
1348      IF(ND.LT.10) THEN
1349       WRITE(DET1(1:1),'(I1.1)')ND
1350       FNAME=MSDATA//DET1
1351      ELSE IF(ND.GE.10.AND.ND.LT.100) THEN
1352       WRITE(DET2(1:2),'(I2.2)')ND
1353       FNAME=MSDATA//DET2
1354      ELSE
1355       WRITE(DET3(1:3),'(I3.3)')ND
1356       FNAME=MSDATA//DET3
1357      END IF
1358     
1359      WRITE(6,2) FNAME
1360    2  FORMAT(' ',' MULTIPLE SCATTERING READ FROM FILE;',80A)
1361      OPEN(UNIT=1,FILE=FNAME,STATUS='OLD')
1362
1363* Read multiple scattering data from file.       
1364      read(1,*) nsmax
1365      write(6,*) ' number of orders of scattering=',nsmax
1366     
1367      i=0
1368   10 i=I+1
1369      read(1,*,end=20) x(i),(d(i,j),j=1,nsmax)
1370      read(1,*,end=20) x(i),(de(i,j),j=1,nsmax)
1371      go to 10
1372   20 CONTINUE
1373      NPTS=i-1
1374      close(1)
1375      WRITE(6,*) ' NPTS=',NPTS,' LPT=',LPT
1376 
1377* Calculate ratio of multiple scattering to total scattering
1378
1379      do j=1,nsmax
1380       sum(j)=0.0
1381       do i=1,npts
1382        sum(j)=sum(j)+d(i,j)
1383        TOTSC(I)=TOTSC(I)+D(I,J)
1384        IF(J.GT.1)MSC(I)=MSC(I)+D(I,J)
1385       end do
1386       write(6,*) ' area of neutrons scattered',j,' times=',sum(j)
1387      end do
1388
1389      SUMT=0.0
1390      SUMM=0.0
1391      SUMD=0.0
1392      do i=1,lpt
1393       XV=XI(i)
1394       WORK1(I)=YV(XV,X,TOTSC,NPTS)
1395       WORK2(I)=YV(XV,X,MSC,NPTS)
1396       SUMT=SUMT+WORK1(I)
1397       SUMM=SUMM+WORK2(I)
1398       SUMD=SUMD+YI(I)
1399      END DO
1400      WRITE(6,*) ' SUMT=',SUMT,' SUMM=',SUMM,' SUMD=',SUMD
1401
1402      DO I=1,LPT
1403       RMS(I)=WORK2(I)*SUMD/SUMT
1404      END DO
1405      END
1406
1407* CALCULATE Y GIVEN X AND ARRAYS OF X,Y VALUES. LINEAR INTERPOLATION.
1408      FUNCTION YV(XV,X,Y,N)
1409      REAL X(N),Y(N)
1410
1411      IF(XV.LT.X(1).OR.XV.GT.X(N)) THEN
1412       YV=0.0
1413       RETURN
1414      END IF
1415
1416      DO I=1,N-1
1417      IF(X(I).LE.XV.AND.XV.LT.X(I+1)) THEN
1418      YV=Y(I)+(XV-X(I))*(Y(I+1)-Y(I))/(X(I+1)-X(I))
1419*      write(6,*) ' xv=',xv,' x(i)=',x(i),' x(i+1)=',x(i+1)
1420*      write(6,*) ' yv=',yv,' y(i)=',y(i),' y(I+1)=',y(I+1)
1421      RETURN
1422      END IF
1423      END DO
1424
1425      END
1426
1427* Error function. Taken from Numerical Recipes p164
1428      FUNCTION ERF(X)
1429      Z=ABS(X)
1430      T=1.0/(1.0+0.5*Z)
1431      ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+
1432     $     T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+
1433     $     T*(1.48851587+T*(-.82215223+T*.17087277)))))))))
1434      IF(X.LT.0) ERFCC=2.-ERFCC
1435      ERF=1.0-ERFCC
1436      RETURN
1437      END
1438 
1439* Calculate time of flight spectrum.
1440* NM=number of masses, M=mass, XS=amplitude, SIGY=width in y.
1441* TRES=resolution fn in tof, L0,L1, lengths in metres, TH angle in rad.
1442* T0 time delay in usec, E1 analyser energy in meV.
1443      SUBROUTINE CTSF(CT,NPTS,NPEAKS,P)
1444      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
1445     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
1446      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
1447      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
1448     $ DTE(20,3),TR0(20),NPEAKS1
1449      REAL CT(NPTS),M
1450      REAL RES(3000),CTM(3000),P(30),AMP(30),WIDTH(30)
1451*       open(unit=4,file='temp.dat',status='new')
1452      NP=0
1453      DO I=1,NPEAKS
1454       IF(XS(I).EQ.0.0) THEN
1455        NP=NP+1
1456        AMP(I)=ABS(P(NP))
1457       ELSE
1458        AMP(I)=ABS(P(NPAR+1)*XS(I))
1459       END IF
1460       IF(WID(I).EQ.0.0)THEN
1461        NP=NP+1
1462        WIDTH(I)=P(NP)
1463       ELSE
1464        WIDTH(I)=WID(I)
1465       END IF
1466      END DO
1467
1468      DO IT=1,NPTS
1469       CT(IT)=0.0
1470      END DO
1471
1472      DO IE=1,1
1473       IF(IE.EQ.1) FACT=1.0
1474       IF(IE.EQ.2) FACT=0.516
1475       IF(IE.EQ.3) FACT=0.079
1476      DO IM=1,NPEAKS
1477*       write(4,*) ' e1=',e1(ie),' M=',am(im)
1478       IF(AM(IM).LT.3) THEN ! Harmonic Oscillator
1479        DSQV=12.0*4.18036*WIDTH(IM)**4/AM(IM)
1480       ELSE ! Debye Approximation
1481        DSQV=12.8*4.18036*WIDTH(IM)**4/AM(IM)
1482       END IF
1483       if(nfse.eq.0) dsqv=0.0
1484       DO IT=1,NPTS
1485        YV=YA(IM,IT,IE)
1486        YV2=YA(IM,IT+1,IE)
1487        if(xi(it).eq.0.0) return
1488        DYDT=(yv2-yv)/(xi(it+1)-xi(it))
1489         if(dydt.eq.0.0) then
1490         YV=YA(IM,IT-1,IE)
1491         YV2=YA(IM,IT,IE)
1492         DYDT=(yv2-yv)/(xi(it+1)-xi(it))
1493         end if
1494        E0=E0A(IT,IE)
1495        Q=QA(IT,IE)
1496        IF(IE.EQ.1) THEN
1497         CPV=CP(YV,Q,AM(IM),WIDTH(IM),DSQV)
1498        ELSE
1499         if(dydt.eq.0.0) then
1500          write(6,*) ' e1=',e1(ie),' m=',am(im)
1501          write(6,*) ' yv2=',yv2,' yv=',yv
1502          write(6,*) ' t2=',xi(it+1),' t=',xi(it)
1503         end if
1504         if(dydt.ne.0.0) then
1505          CPV=CP2(YV,YV2,Q,AM(IM),WIDTH(IM),DSQV)/dydt
1506         else
1507          CPV=0.0
1508         end if
1509         if(it.eq.npts) cpv=0.0
1510*         if(xi(it).eq.110.and.IE.eq.2) then
1511*          CPV=CP2W(YV,YV2,Q,AM(IM),WIDTH(IM),DSQV)/dydt
1512*         end if
1513        END IF
1514        CTM(IT)=AM(IM)*AMP(IM)*CPV*E0**0.1/Q
1515*        write(4,11) xi(it),yv,cpv,ctm(it)
1516   11    Format(' ',' t=',1pe11.4,' yv=',1pe11.4,' cpv=',1pe11.4,
1517     $' ctm=',1pe11.4)       
1518        RES(IT)=TRES(IM,IT,IE)
1519       END DO
1520       CALL C06EKE(1,CTM,RES,NPTS,IFAIL) ! Convolve with resn fn.
1521       DO IT=1,NPTS
1522        CT(IT)=CT(IT)+FACT*CTM(IT)
1523       END DO
1524
1525      END DO
1526      END DO
1527
1528      IF(NBACK.EQ.1) THEN
1529       N=NP+4
1530       DO IT=1,NPTS
1531        X=XI(IT)
1532       CT(IT)=CT(IT)+P(N-2)*(TMAX-X)+P(N-1)*(TMAX-X)**2
1533     $ +P(N)*(TMAX-X)**3
1534       END DO                   
1535      END IF
1536
1537      IF(NMSSUB.EQ.1.AND.NFACT.EQ.1) THEN ! Muliple scattering correction.
1538       DO I=1,NPTS
1539        CT(I)=CT(I)+P(N+1)*RMS(I)
1540       END DO
1541      ELSE IF(NMSSUB.EQ.1.AND.NFACT.NE.1) THEN
1542      DO I=1,NPTS
1543       CT(I)=CT(I)+RMS(I)
1544      END DO
1545      END IF
1546
1547      close(4)
1548      END
1549
1550* Calculates J(y) including first Sears Correction terms
1551      FUNCTION CP2W(Y1,y2,Q,M,SIG1,D2V)
1552      REAL Y,M,J0,J3
1553
1554      NP=101 ! Number of points in integration
1555      CP=0.0
1556      SIG=ABS(SIG1)
1557      DY=(y2-y1)/(np-1)
1558*       write(4,*) ' y1=',y1,' y2=',y2,' dy=',dy
1559      do i=1,np
1560       y=y1+dy*(i-1)
1561       X=Y/SIG
1562       J0=EXP(-X**2/2)/(2.506628*SIG) ! 2.5066=SQRT(2PI)
1563       D3J=X*(3-X**2)*J0/SIG**3
1564       J3=-M*D2V*D3J/(150.49*Q)
1565       if(d2v.ne.0.0) then
1566        CP=CP+J0+J3
1567       else
1568        cp=cp+j0
1569       end if
1570*       write(4,*) ' y=',y,' j0=',j0,' cp=',cp
1571      END DO
1572      CP2W=CP*dy
1573
1574      END
1575
1576* SUBROUTINE TO PRODUCE OUTPUT FILE CONTAINING FITTED
1577* PARAMETERS
1578      SUBROUTINE FOUT(NM,NDETMIN,ndet,
1579     $AM,WT,WTE,SIG,SIGE,CHI,FIN,ipno)
1580      PARAMETER (NP=20,nd=200)
1581      REAL wt(np,nd),wte(np,nd),sig(np,nd)
1582     *,sige(np,nd),Chi(nd),am(nm)
1583      CHARACTER FIN*70,fileout*80
1584     
1585*      write(6,*) ' nm=',nm,' ndetmin=',ndetmin,' ndet=',ndet
1586      ndetmax=ndetmin+ndet-1
1587     
1588*      do i=ndetmin,ndetmax
1589*      write(6,*) i, wt(1,i),wte(1,i)
1590*      end do
1591     
1592     
1593      FileOUT=FIN//'P.dat'
1594      OPEN(UNIT=4,FILE=fileout,STATUS='NEW')
1595      WRITE(4,*) NM
1596      NDTOT=NDETMAX-NDETMIN+1
1597      WRITE(4,*) NDTOT
1598      WRITE(4,*) ' '
1599      DO IM=1,NM
1600       WRITE(4,*) ' '
1601       WRITE(4,*) am(IM)
1602       WRITE(4,*) ' '
1603      WRITE(4,*) '   ','      WEIGHT    ','     ERROR     ',
1604     $'     SIGMA     ','     ERROR     ','   CHI-SQ      '
1605       WRITE(6,*) ' '
1606       WRITE(6,*) ' FOR MASS', am(IM)
1607       WRITE(6,*) ' '
1608      WRITE(6,*) '   ','      WEIGHT    ','     ERROR     ',
1609     $'     SIGMA     ','     ERROR     ','   CHI-SQ      '
1610       DO ID=Ndetmin,Ndetmax
1611        WRITE(6,15) ID,WT(Im,Id),WTE(Im,Id),SIG(Im,Id)
1612     $  ,SIGE(Im,Id),CHI(ID)
1613        WRITE(4,15) ID,WT(Im,Id),WTE(Im,Id),SIG(Im,Id)
1614     $  ,SIGE(Im,Id),CHI(ID)
1615   15   FORMAT(' ',I3,1P5E15.4)
1616       END DO
1617      END DO
1618      write(4,*) ipno
1619      CLOSE(4)
1620
1621      end
1622       
1623
1624* Calculate time of flight spectrum.
1625* NM=number of masses, M=mass, XS=amplitude, SIGY=width in y.
1626* TRES=resolution fn in tof, L0,L1, lengths in metres, TH angle in rad.
1627* T0 time delay in usec, E1 analyser energy in meV.
1628      SUBROUTINE CTSF2(CT2,NPTS,NPEAKS,P,N)
1629      COMMON QA(3000,3),E0A(3000,3),ya(20,3000,3),TMAX
1630     $,TRES(20,3000,3),E1(3),NBACK,NMSSUB,NFACT,NFSE
1631      COMMON AM(20),DTG(20,3),XS(20),WID(20),NPAR
1632      COMMON XI(3000),YI(3000),EI(3000),V(30),RMS(3000),
1633     $ DTE(20,3),TR0(20),NPEAKS1
1634      REAL CT2(20,1024),M,ct(npts),amp1(30)
1635      REAL RES(3000),CTM(3000),P(30),AMP(30),WIDTH(30)
1636
1637      do ip=1,npeaks ! Start of loop over peaks
1638      NP=0
1639      DO I=1,NPEAKS
1640       IF(XS(I).EQ.0.0) THEN
1641        NP=NP+1
1642        AMP1(I)=ABS(P(NP))
1643*        write(6,*) i,' xs=0.amp1=',amp1(i)
1644       ELSE
1645        AMP1(I)=ABS(P(NPAR+1)*XS(I))
1646*        write(6,*) i,' xs.ne.0.amp1=',amp1(i)
1647       END IF
1648       
1649       IF(WID(I).EQ.0.0)THEN
1650        NP=NP+1
1651        WIDTH(I)=P(NP)
1652       ELSE
1653        WIDTH(I)=WID(I)
1654       END IF
1655*       write(6,*) i, ' amp1=',amp1(i),' width=',width(i)
1656      END DO
1657
1658      do i=1,npeaks
1659      if(i.ne.ip) amp(i)=0.0
1660      if(i.eq.ip) amp(i)=amp1(i)
1661      end do
1662*      write(6,*) ip, amp(1),amp(2),amp(3)
1663
1664      DO IT=1,NPTS
1665       CT(IT)=0.0
1666      END DO
1667
1668      DO IE=1,1
1669       IF(IE.EQ.1) FACT=1.0
1670       IF(IE.EQ.2) FACT=0.516
1671       IF(IE.EQ.3) FACT=0.079
1672      DO IM=1,NPEAKS
1673*       write(6,*) ' e1=',e1(ie),' M=',am(im)
1674       IF(AM(IM).LT.3) THEN ! Harmonic Oscillator
1675        DSQV=12.0*4.18036*WIDTH(IM)**4/AM(IM)
1676       ELSE ! Debye Approximation
1677        DSQV=12.8*4.18036*WIDTH(IM)**4/AM(IM)
1678       END IF
1679       if(nfse.eq.0) dsqv=0.0
1680       DO IT=1,NPTS
1681        YV=YA(IM,IT,IE)
1682        YV2=YA(IM,IT+1,IE)
1683        if(xi(it).eq.0.0) return
1684        DYDT=(yv2-yv)/(xi(it+1)-xi(it))
1685         if(dydt.eq.0.0) then
1686         YV=YA(IM,IT-1,IE)
1687         YV2=YA(IM,IT,IE)
1688         DYDT=(yv2-yv)/(xi(it+1)-xi(it))
1689         end if
1690        E0=E0A(IT,IE)
1691        Q=QA(IT,IE)
1692        IF(IE.EQ.1) THEN
1693         CPV=CP(YV,Q,AM(IM),WIDTH(IM),DSQV)
1694        ELSE
1695         if(dydt.eq.0.0) then
1696          write(6,*) ' e1=',e1(ie),' m=',am(im)
1697          write(6,*) ' yv2=',yv2,' yv=',yv
1698          write(6,*) ' t2=',xi(it+1),' t=',xi(it)
1699         end if
1700         if(dydt.ne.0.0) then
1701          CPV=CP2(YV,YV2,Q,AM(IM),WIDTH(IM),DSQV)/dydt
1702         else
1703          CPV=0.0
1704         end if
1705         if(it.eq.npts) cpv=0.0
1706        END IF
1707        CTM(IT)=AM(IM)*AMP(IM)*CPV*E0**0.1/Q
1708*        if(im.eq.3.and.ip.eq.3) write(6,11) xi(it),yv,cpv,ctm(it)
1709   11    Format(' ',' t=',1pe11.4,' yv=',1pe11.4,' cpv=',1pe11.4,
1710     $' ctm=',1pe11.4)       
1711        RES(IT)=TRES(IM,IT,IE)
1712       END DO
1713       CALL C06EKE(1,CTM,RES,NPTS,IFAIL) ! Convolve with resn fn.
1714       DO IT=1,NPTS
1715        CT(IT)=CT(IT)+FACT*CTM(IT)
1716       END DO
1717
1718      END DO
1719      END DO
1720
1721      IF(NBACK.EQ.1) THEN
1722       N=NP+4
1723       DO IT=1,NPTS
1724        X=XI(IT)
1725       CT(IT)=CT(IT)+P(N-2)*(TMAX-X)+P(N-1)*(TMAX-X)**2
1726     $ +P(N)*(TMAX-X)**3
1727       END DO                   
1728      END IF
1729
1730      IF(NMSSUB.EQ.1.AND.NFACT.EQ.1) THEN ! Muliple scattering correction.
1731       DO I=1,NPTS
1732        CT(I)=CT(I)+P(N+1)*RMS(I)
1733       END DO
1734      ELSE IF(NMSSUB.EQ.1.AND.NFACT.NE.1) THEN
1735      DO I=1,NPTS
1736       CT(I)=CT(I)+RMS(I)
1737      END DO
1738      END IF
1739
1740      do i=1,npts
1741      ct2(ip,i)=ct(i)
1742      end do
1743     
1744      end do ! end of loop over peaks
1745
1746
1747      END