Ticket #8738: QLwat_subs.f90

File QLwat_subs.f90, 18.1 KB (added by Samuel Jackson, 7 years ago)

Fortran code for QLWater

Line 
1      SUBROUTINE init_paras
2      INCLUDE 'res_par.f90'
3      COMMON /FFTCOM/ FRES(m_d2),FWRK(m_d2),XJ(m_d),TWOPIK(m_d1),NFFT
4      COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT
5      COMMON /DINTRP/ IPDAT(m_d),XPDAT(m_d)
6      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
7      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
8      COMMON /GRDCOM/ DDDPAR(m_d,m_p),FR2PIK(m_d2,2)
9      COMMON /QW1COM/ QW1(m_sp),SIGQW1(m_sp),ISPEC
10      COMMON /WRKCOM/ WORK(m_d2,2)
11      do n=1,m_d
12       XJ(n)=0.0
13       XDAT(n)=0.0
14       DAT(n)=0.0
15       SIG(n)=0.0
16       IPDAT(n)=0
17       XPDAT(n)=0.0
18       FIT(n)=0.0
19       RESID(n)=0.0
20       do m=1,m_p
21        DDDPAR(n,m)=0.0
22       end do
23      end do
24      do n=1,m_d1
25       TWOPIK(n)=0.0
26       do m=1,3
27        EXPF(n,m)=0.0
28       end do
29      end do
30      do n=1,m_d2
31       FRES(n)=0.0
32       FWRK(n)=0.0
33       do m=1,2
34        FR2PIK(n,m)=0.0
35        WORK(n,m)=0.0
36       end do
37      end do
38      do n=1,m_p
39       FITP(n)=0.0
40       do m=1,2
41        SCLVEC(n,m)=0.0
42       end do
43      end do
44      do n=1,m_sp
45       QW1(n)=0.0
46       SIGQW1(n)=0.0
47      end do
48      END
49C     ------------------------------------
50      SUBROUTINE OUTPRM(P,C,NP,NFEW,CNORM)
51      INCLUDE 'mod_files.f90'
52      REAL P(*),C(NP,*)
53      IF (NFEW.NE.1) RETURN
54       OPEN(UNIT=1,FILE=fileout1,STATUS='old',FORM='formatted',
55     1 access='append')
56      WRITE(NFEW,100) P(3),P(1),P(2),P(4)
57      WRITE(NFEW,100) P(5),P(6),P(7)
58      CSCALE=2.0*CNORM
59      do J=1,NP
60       do I=1,NP
61        C(I,J)=CSCALE*C(I,J)
62       end do
63      end do
64      WRITE(NFEW,100) C(3,3)
65      WRITE(NFEW,100) C(3,5),C(5,5)
66      WRITE(NFEW,100) C(3,6),C(5,6),C(6,6)
67      WRITE(NFEW,100) C(3,7),C(5,7),C(6,7),C(7,7)
68 100  FORMAT(1PE13.4,4E13.4)
69      WRITE(NFEW,*)' -------------------------------------------------'
70      close(unit=1)
71      return
72      END
73c
74C***<calculate data & chi-squared>**************************************
75      FUNCTION CCHI(V)
76      INCLUDE 'res_par.f90'
77      INCLUDE 'options.f90'
78      COMMON /FFTCOM/ FRES(m_d2),FWRK(m_d2),XJ(m_d),TWOPIK(m_d1),NFFT
79      COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT
80      COMMON /DINTRP/ IPDAT(m_d),XPDAT(m_d)
81      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
82      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
83      COMMON /GRDCOM/ DDDPAR(m_d,m_p),FR2PIK(m_d2,2)
84      COMMON /QW1COM/ QW1(m_sp),SIGQW1(m_sp),ISPEC
85      COMMON /WRKCOM/ WORK(m_d2,2)
86      COMMON /BESSJ/  AJ0,AJ1,AJ2
87      REAL   V(*)
88      DATA   E3,E9 /3.0,9.0/
89      PI=3.141592654
90      VSMALL=0.0001
91      small=0.001
92      CCHI=1.0E+10
93      CHI=0.0
94      BNRM=1.0
95      B1=BSCL*V(1)
96      B2=BSCL*V(2)
97      A0=ASCL*V(3)
98      DELTAX=V(4)
99      NFT2=NFFT/2+1
100      CALL CXSHFT(FRES,DELTAX,TWOPIK,FR2PIK,FR2PIK(1,2),NFT2)
101      CALL VRFILL(WORK,A0,NFT2)
102      XNSCL=-LOG(1.0E-7)/(TWOPIK(2)-TWOPIK(1))
103      IF (NFEW.EQ.1) THEN
104        E3B=1.0/E3
105        E9B=1.0/E9
106        A1=ASCL*V(5)
107        SIG1=WSCL*V(6)/GSCL
108        SIG2=WSCL*V(7)/(GSCL*3.0)
109        IF (SIG1.LT.0.0 .OR. SIG2.LT.0.0) RETURN
110        CALL VRFILL(EXPF(1,1),0.0,NFT2)
111        CALL VRFILL(EXPF(1,2),0.0,NFT2)
112        CALL VRFILL(EXPF(1,3),0.0,NFT2)
113        NXFT21=1+NINT(XNSCL/(ABS(SIG1)+1.0E-10))
114        NXFT22=1+NINT(XNSCL/(ABS(SIG2)+1.0E-10))
115        NXFT2=MIN(MAX(NXFT21,NXFT22),NFT2)
116        do I=1,NXFT2
117          EXPIJ1=EXP(-TWOPIK(I)*SIG1)
118          EXPIJ2=EXP(-TWOPIK(I)*SIG2)
119          EXPIJ=EXPIJ1*(AJ0+EXPIJ2*(AJ1+AJ2*EXPIJ2**2))
120          WORK(I,1)=WORK(I,1)+A1*EXPIJ
121          EXPF(I,1)=EXPIJ
122          EXPF(I,2)=E3B*EXPIJ1*EXPIJ2*(AJ1+E3*AJ2*EXPIJ2**2)
123          EXPF(I,3)=E9B*EXPIJ1*EXPIJ2*(AJ1+E9*AJ2*EXPIJ2**2)
124        end do
125      ENDIF
126      CALL VMLTRC(WORK,FR2PIK,NFT2,FWRK)
127      CALL FOUR2(FWRK,NFFT,1,-1,-1)
128      CALL DEGRID(FWRK,FIT)
129      X1=XDAT(1)
130      if(o_bgd.eq.2) BNRM=(B2-B1)/(XDAT(NDAT)-X1)
131C       avoid conflict BNORM with ModPars
132      do I=1,NDAT
133        FIT(I)=FIT(I)+B1
134         if(o_bgd.eq.2) FIT(I)=FIT(I)+BNRM*(XDAT(I)-X1)
135        DIF=FIT(I)-DAT(I)
136        RESID(I)=DIF*SIG(I)
137        CHI=CHI+DIF*RESID(I)
138      end do
139      CCHI=CHI/(2.0*FLOAT(NDAT))
140      sm2=small*small
141      v12=V(1)-V(2)
142      if(v12.lt.vsmall)v12=1.0
143      if(o_bgd.le.1) CCHI=CCHI+v12*v12/sm2
144      RETURN
145      END
146C
147C**<search for one more & refine amplitudes>****************************
148C
149      SUBROUTINE SEARCHw(GRAD,HESS,DPAR,NFEW,INDX,COVAR,FITP)
150      INCLUDE 'res_par.f90'
151      INCLUDE 'options.f90'
152      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
153      COMMON /QW1COM/ QW1(m_sp),SIGQW1(m_sp),ISPEC
154      REAL     GRAD(*),HESS(*),DPAR(*),COVAR(*),FITP(*)
155      INTEGER  INDX(*)
156      J=4+3*NFEW
157      DXLOG=0.85
158      NSRCH=NINT(LOG(5.0*GSCL/WSCL)/LOG(DXLOG))
159      CMIN=1.0E20
160      SIG1=0.
161      SIG2=0.
162      FITP(J-1)=1.0
163      do I2=1,NSRCH
164        FITP(J)=1.0
165        FITP(J-2)=0.5
166        do I1=1,NSRCH
167          CALL REFINA(GRAD,HESS,DPAR,3+NFEW,DETLOG,INDX,COVAR)
168          CNORM=CCHI(FITP)
169          IF (CNORM.LT.CMIN) THEN
170            CMIN=CNORM
171            SIG1=FITP(J-1)
172            SIG2=FITP(J)
173          ENDIF
174          FITP(J)=FITP(J)*DXLOG
175        end do
176        FITP(J-1)=FITP(J-1)*DXLOG
177      end do
178      FITP(J-1)=SIG1
179      FITP(J)=SIG2
180      CALL REFINA(GRAD,HESS,DPAR,3+NFEW,DETLOG,INDX,COVAR)
181      END
182C**<refinement>*********************************************************
183      SUBROUTINE REFINAw(GRAD,HESS,DPAR,NP,DETLOG,INDX,COVAR)
184      INCLUDE 'res_par.f90'
185      INCLUDE 'options.f90'
186      COMMON /FFTCOM/ FRES(m_d2),FWRK(m_d2),XJ(m_d),TWOPIK(m_d1),NFFT
187      COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT
188      COMMON /DINTRP/ IPDAT(m_d),XPDAT(m_d)
189      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
190      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
191      COMMON /GRDCOM/ DDDPAR(m_d,m_p),FR2PIK(m_d2,2)
192      REAL            GRAD(*),HESS(NP,*),DPAR(*),COVAR(NP,*)
193      INTEGER         INDX(*)
194      NFT2=NFFT/2+1
195      CNORM=CCHI(FITP)
196      CALL VRFILL(HESS,0.0,NP*NP)
197      CALL VCOPY(FR2PIK,FWRK,NFFT+2)
198      CALL FOUR2(FWRK,NFFT,1,-1,-1)
199      CALL DEGRID(FWRK,DDDPAR(1,3))
200      if (NFEW.eq.1)then
201        CALL VMLTRC(EXPF(1,1),FR2PIK,NFT2,FWRK)
202        CALL FOUR2(FWRK,NFFT,1,-1,-1)
203        CALL DEGRID(FWRK,DDDPAR(1,4))
204      endif
205      CALL GRADPR(GRAD,RESID,NDAT,NP,SCLVEC)
206      CALL HESS1(HESS,NP,SCLVEC,0.3,NFEW)
207      CALL INVERT(HESS,COVAR,NP,INDX,DETLOG)
208      CALL NEWEST(COVAR,GRAD,NP,NFEW,DPAR,FITP)
209      CNORM=CCHI(FITP)
210      CALL GRADPR(GRAD,RESID,NDAT,NP,SCLVEC)
211      CALL NEWEST(COVAR,GRAD,NP,NFEW,DPAR,FITP)
212      END
213C     -------------------------------------------------------
214      SUBROUTINE REFINEw(GRAD,HESS,NP,DETLOG,INDX,COVAR,STEPSZ)
215      INCLUDE 'res_par.f90'
216      INCLUDE 'options.f90'
217      COMMON /FFTCOM/ FRES(m_d2),FWRK(m_d2),XJ(m_d),TWOPIK(m_d1),NFFT
218      COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT
219      COMMON /DINTRP/ IPDAT(m_d),XPDAT(m_d)
220      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
221      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
222      COMMON /WRKCOM/ WORK(m_d2,2)
223      COMMON /GRDCOM/ DDDPAR(m_d,m_p),FR2PIK(m_d2,2)
224      COMMON /QW1COM/ QW1(m_sp),SIGQW1(m_sp),ISPEC
225      COMMON /BESSJ/  AJ0,AJ1,AJ2
226      REAL GRAD(*),HESS(NP,*),COVAR(NP,*)
227      INTEGER INDX(*)
228C
229      NFT2=NFFT/2+1
230      CNORM=CCHI(FITP)
231      CALL VRFILL(HESS,0.0,NP*NP)
232      CALL VMLTRC(WORK,FR2PIK(1,2),NFT2,FWRK)
233      CALL VMLTRC(TWOPIK,FWRK,NFT2,WORK)
234      CALL VMLTIC(WORK,NFT2,WORK)
235      CALL FOUR2(FWRK,NFFT,1,-1,-1)
236      CALL DEGRID(FWRK,DDDPAR(1,4))
237      CALL FOUR2(WORK,NFFT,1,-1,-1)
238      CALL DEGRID(WORK,FWRK)
239      CALL VRDOTR(RESID,FWRK,NDAT,HESS(4,4))
240      CALL VCOPY(FR2PIK,FWRK,NFFT+2)
241      CALL FOUR2(FWRK,NFFT,1,-1,-1)
242      CALL DEGRID(FWRK,DDDPAR(1,3))
243      CALL VCOPY(FR2PIK(1,2),FWRK,NFFT+2)
244      CALL FOUR2(FWRK,NFFT,1,-1,-1)
245      CALL DEGRID(FWRK,WORK)
246      CALL VRDOTR(RESID,WORK,NDAT,HESS(3,4))
247      HESS(4,3)=HESS(3,4)
248      IF (NFEW.EQ.1) THEN
249        A1=FITP(5)*ASCL
250        CALL VMLTRC(EXPF(1,1),FR2PIK,NFT2,WORK)
251        CALL VCOPY(WORK,FWRK,NFFT+2)
252        CALL FOUR2(FWRK,NFFT,1,-1,-1)
253        CALL DEGRID(FWRK,DDDPAR(1,5))
254        CALL VMLTRC(TWOPIK,WORK,NFT2,FWRK)
255        CALL VCOPY(FWRK,WORK,NFFT+2)
256        CALL FOUR2(FWRK,NFFT,1,-1,-1)
257        CALL DEGRID(FWRK,DDDPAR(1,6))
258        CALL HESS0(HESS,NP,RESID,NDAT,DDDPAR(1,6),A1,5)
259        CALL VMLTIC(WORK,NFT2,WORK)
260        CALL VCOPY(WORK,FWRK,NFFT+2)
261        CALL FOUR2(FWRK,NFFT,1,-1,-1)
262        CALL DEGRID(FWRK,WORK(1,2))
263        CALL VRDOTR(RESID,WORK(1,2),NDAT,HESS(4,5))
264        HESS(5,4)=HESS(4,5)
265        CALL VMLTRC(TWOPIK,WORK,NFT2,FWRK)
266        CALL VCOPY(FWRK,WORK,NFFT+2)
267        CALL FOUR2(FWRK,NFFT,1,-1,-1)
268        CALL DEGRID(FWRK,WORK(1,2))
269        CALL VRDOTR(RESID,WORK(1,2),NDAT,SM)
270        HESS(4,6)=-A1*SM
271        HESS(6,4)=HESS(4,6)
272        CALL VMLTIC(WORK,NFT2,WORK)
273        CALL FOUR2(WORK,NFFT,1,-1,-1)
274        CALL DEGRID(WORK,FWRK)
275        CALL VRDOTR(RESID,FWRK,NDAT,SM)
276        HESS(6,6)=-A1*SM
277        CALL VMLTRC(EXPF(1,2),FR2PIK,NFT2,WORK)
278        CALL VMLTRC(TWOPIK,WORK,NFT2,FWRK)
279        CALL VCOPY(FWRK,WORK,NFFT+2)
280        CALL FOUR2(FWRK,NFFT,1,-1,-1)
281        CALL DEGRID(FWRK,DDDPAR(1,7))
282        CALL HESS0(HESS,NP,RESID,NDAT,DDDPAR(1,7),A1,5)
283        CALL VMLTIC(WORK,NFT2,WORK)
284        CALL VMLTRC(TWOPIK,WORK,NFT2,FWRK)
285        CALL VCOPY(FWRK,WORK,NFFT+2)
286        CALL FOUR2(FWRK,NFFT,1,-1,-1)
287        CALL DEGRID(FWRK,WORK(1,2))
288        CALL VRDOTR(RESID,WORK(1,2),NDAT,SM)
289        HESS(4,7)=-A1*SM
290        HESS(7,4)=HESS(4,7)
291        CALL VMLTIC(WORK,NFT2,WORK)
292        CALL FOUR2(WORK,NFFT,1,-1,-1)
293        CALL DEGRID(WORK,FWRK)
294        CALL VRDOTR(RESID,FWRK,NDAT,SM)
295        HESS(6,7)=-A1*SM
296        HESS(7,6)=HESS(6,7)
297        CALL VMLTRC(EXPF(1,3),FR2PIK,NFT2,WORK)
298        CALL VMLTRC(TWOPIK,WORK,NFT2,FWRK)
299        CALL VMLTRC(TWOPIK,FWRK,NFT2,WORK)
300        CALL FOUR2(WORK,NFFT,1,-1,-1)
301        CALL DEGRID(WORK,FWRK)
302        CALL VRDOTR(RESID,FWRK,NDAT,SM)
303        HESS(7,7)=A1*SM
304      ENDIF
305      CALL GRADPR(GRAD,RESID,NDAT,NP,SCLVEC(1,2))
306      CALL HESS1(HESS,NP,SCLVEC(1,2),STEPSZ,0)
307      CALL INVERT(HESS,COVAR,NP,INDX,DETLOG)
308      RETURN
309      END
310C***<see the fit>*******************************************************
311      SUBROUTINE SEEFITw(SIGPAR,CNORM,PRMSV,SIGSV)
312      INCLUDE 'res_par.f90'
313      COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT
314      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
315      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
316      COMMON /WRKCOM/ WORK(m_d2,2)
317      INCLUDE 'mod_files.f90'
318      REAL       SIGPAR(*),PRMSV(*),SIGSV(*)
319      OPEN(UNIT=53,FILE=lptfile,STATUS='old',FORM='formatted',
320     1 access='append')
321      ERRSCL=SQRT(CNORM)
322      NQUASI=0
323      IF (NFEW.EQ.1) NQUASI=3
324      WRITE(53,100) NQUASI
325 100  FORMAT(' Best-fit assuming no. of quasi-elastic lines = ',I2)
326      WRITE(53,120) CNORM
327 120  FORMAT(' Normalised Chi-squared = ',F11.4)
328      WRITE(53,130) FITP(1)*BSCL,SIGPAR(1)*BSCL*ERRSCL
329 130  FORMAT(' Background(Xmin) = ',1PE12.3,'  +- ',E10.2)
330      WRITE(53,99) SIGPAR(1),BSCL,ERRSCL
331  99  FORMAT(' Sigpar = ',3(1PE12.3,1x))
332      WRITE(53,140) FITP(2)*BSCL,SIGPAR(2)*BSCL*ERRSCL
333 140  FORMAT(' Background(Xmax) = ',1PE12.3,'  +- ',E10.2)
334      WRITE(53,150) FITP(4)*GSCL*1000.0,
335     1  SIGPAR(4)*GSCL*ERRSCL*1000.0
336 150  FORMAT(' Zero offset       = ',F12.2,'  +- ',F10.2,'  ueV')
337      WRITE(53,*)' Elastic line'
338      WRITE(53,160) FITP(3)*ASCL,SIGPAR(3)*ASCL*ERRSCL
339 160  FORMAT(5X,' Amplitude  =   ',1PE13.4,'  +- ',E11.3)
340      PRMSV(1)=FITP(3)*ASCL
341      SIGSV(1)=SIGPAR(3)*ASCL*ERRSCL
342      IF (NFEW.EQ.1) THEN
343        WRITE(53,161)
344 161    FORMAT(' Quasi-elastic parameters')
345        WRITE(53,170) 2.0*FITP(6)*WSCL,2.0*SIGPAR(6)*WSCL*ERRSCL
346 170    FORMAT(5X,' FWHM0      =   ',F13.2,'  +- ',F11.2)
347        WRITE(53,171) 2.0*FITP(7)*WSCL,2.0*SIGPAR(7)*WSCL*ERRSCL
348 171    FORMAT(5X,' FWHM1      =   ',F13.2,'  +- ',F11.2)
349        WRITE(53,180) FITP(5)*ASCL,SIGPAR(5)*ASCL*ERRSCL
350 180    FORMAT(5X,' Amplitude  =   ',1PE13.4,'  +- ',E11.3)
351        PRMSV(2)=FITP(5)*ASCL
352        SIGSV(2)=SIGPAR(5)*ASCL*ERRSCL
353        PRMSV(3)=2.0*FITP(6)*WSCL
354        SIGSV(3)=2.0*SIGPAR(6)*WSCL*ERRSCL
355        PRMSV(4)=2.0*FITP(7)*WSCL
356        SIGSV(4)=2.0*SIGPAR(7)*WSCL*ERRSCL
357      ENDIF
358      close(unit=53)
359      END
360C     ----------------------------------
361      SUBROUTINE ERRBARw(COVAR,NP,SIGPAR,COV12)
362      REAL COVAR(NP,*),SIGPAR(*)
363      SMALL=1.0E-20
364      do I=1,NP
365        SIGPAR(I)=SQRT(2.0*ABS(COVAR(I,I))+SMALL)
366      end do
367      IF (NP.EQ.7) COV12=2.0*COVAR(6,7)
368      END
369C     ----------------------------------------------------------
370      SUBROUTINE INTDSPw(X,Y,XL,XH,YL,YH,SIGPAR,ERRSCL,XMIN,XMAX,COV12)
371      INCLUDE 'res_par.f90'
372      COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6)
373      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
374      COMMON /BESSJ/  AJ0,AJ1,AJ2
375      REAL   X(*),Y(*),XL(*),XH(*),YL(*),YH(*),SIGPAR(*)
376      XMAX=0.0
377      CALL VRFILL(X,0.0,2*(NFEW+1))
378      CALL VRFILL(Y,0.0,2*(NFEW+1))
379      Y(2)=FITP(3)*ASCL
380      XL(1)=0.0
381      XH(1)=0.0
382      YL(1)=Y(2)-ASCL*SIGPAR(3)*ERRSCL
383      YH(1)=Y(2)+ASCL*SIGPAR(3)*ERRSCL
384      IF (NFEW.EQ.1) THEN
385       A1=FITP(5)*ASCL
386       FWHM1=FITP(6)*WSCL*2.0
387       FWHM2=FITP(7)*WSCL*2.0
388       SIGA1=SIGPAR(5)*ASCL*ERRSCL
389       SIGFW1=SIGPAR(6)*WSCL*ERRSCL*2.0
390       SIGFW2=SIGPAR(7)*WSCL*ERRSCL*2.0
391       COV12B=COV12*(WSCL*ERRSCL*2.0)**2
392       CALL VRFILL(X(3),FWHM1,2)
393       Y(4)=A1*AJ0
394       SIGX=SIGFW1
395       SIGY=SIGA1*AJ0
396       XL(2)=X(4)-SIGX
397       XH(2)=X(4)+SIGX
398       YL(2)=Y(4)-SIGY
399       YH(2)=Y(4)+SIGY
400       CALL VRFILL(X(5),FWHM1+FWHM2/3.0,2)
401       Y(6)=A1*AJ1
402       SIGX=SQRT(SIGFW1**2+(SIGFW2**2)/9.0+2.0*COV12B/3.0)
403       SIGY=SIGA1*AJ1
404       XL(3)=X(6)-SIGX
405       XH(3)=X(6)+SIGX
406       YL(3)=Y(6)-SIGY
407       YH(3)=Y(6)+SIGY
408       CALL VRFILL(X(7),FWHM1+FWHM2,2)
409       Y(8)=A1*AJ2
410       SIGX=SQRT(SIGFW1**2+SIGFW2**2+2.0*COV12B)
411       SIGY=SIGA1*AJ2
412       XL(4)=X(8)-SIGX
413       XH(4)=X(8)+SIGX
414       YL(4)=Y(8)-SIGY
415       YH(4)=Y(8)+SIGY
416       XMAX=MAX(XMAX,XH(2),XH(3),XH(4))
417       XMIN=-XMAX/50.0
418       XMAX=XMAX*1.03
419      ELSE
420       XMAX=WSCL/20.0
421       XMIN=-XMAX/5.0
422      ENDIF
423      END
424C     ------------------------------------------------------------------
425      SUBROUTINE PROBNw(CNORM,NDAT,DETLOG,NFEW,PROBLG)
426      INCLUDE 'res_par.f90'
427      INCLUDE 'mod_files.f90'
428      INCLUDE 'options.f90'
429      COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL
430      EXTERNAL FCTNLG
431      CHISQ=CNORM*FLOAT(NDAT)
432      DETLOG=DETLOG-FLOAT(NFEW*2)*LOG10(ASCL*WSCL)
433      CHI2LG=-LOG10(2.7182818)*CHISQ/2.0
434      PROBLG=CHI2LG-(0.5*DETLOG)-(FLOAT(NFEW)*LOG10(ASCL*WSCL))
435      PROBLG=PROBLG+FCTNLG(NFEW)
436      PROBLG=PROBLG+(FLOAT(NFEW)*LOG10(4.0*3.141592654))
437      OPEN(UNIT=53,FILE=lptfile,STATUS='old',FORM='formatted',
438     1 access='append')
439      WRITE(53,110) PROBLG
440 110  FORMAT(' Log10[Prob(Stretched exp|{Data})] = ',F11.1)
441      IF (NFEW.LT.1) WRITE(53,*)' -------------------------'
442      close(unit=53)
443      END
444C     ---------------------------------
445      SUBROUTINE PRBOUTw(P,NP,NQ,POUT)
446      INCLUDE 'res_par.f90'
447      INCLUDE 'mod_files.f90'
448      REAL  P(4,m_sp),POUT(4,m_sp)
449      J=NQ
450      SM=0.0
451      do I=1,NP
452        SM=SM+10.0**(P(I,J)-P(NP,J))
453      end do
454      PLNORM=LOG10(SM)+P(NP,J)
455      do I=1,NP
456        P(I,J)=P(I,J)-PLNORM
457      end do
458      do I=1,NP
459        POUT(I,J)=P(I,J)
460       end do
461      END
462C     ---------------------------------
463      SUBROUTINE BESINT(AJ0,AJ1,AJ2,QA)
464      INCLUDE 'mod_files.f90'
465      BJ0=BESSJ0(QA)
466      BJ1=BESSJ1(QA)
467      BJ2=BJ1*2.0/QA-BJ0
468      AJ0=BJ0**2
469      AJ1=3.0*BJ1**2
470      AJ2=5.0*BJ2**2
471      RETURN
472      END
473C     ---------------------------------
474      SUBROUTINE BTINIT(X,N,XMIN,XMAX,NBMIN)
475      REAL    X(*)
476      INTEGER NBMIN(*)
477      NBMIN(1)=N
478      NBMIN(2)=NINT(0.1*FLOAT(N))
479      DX=(XMAX-XMIN)/FLOAT(N-1)
480      X(1)=XMIN
481      do I=2,N
482       X(I)=X(I-1)+DX
483      end do
484      RETURN
485      END
486C***<Numerical Recipes>*************************************************
487      FUNCTION bessj0(x)
488      REAL bessj0,x
489      REAL ax,xx,z
490      DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
491     *                 s1,s2,s3,s4,s5,s6,y
492      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
493     *     s5,s6
494      DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
495     * -.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
496     * .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
497      DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
498     * 651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
499     * s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
500     * 59272.64853d0,267.8532712d0,1.d0/
501      if(abs(x).lt.8.)then
502        y=x**2
503        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
504     *         (s4+y*(s5+y*s6)))))
505      else
506        ax=abs(x)
507        z=8./ax
508        y=z**2
509        xx=ax-.785398164
510        bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
511     *              p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
512      endif
513      END
514C     ------------------
515      FUNCTION bessj1(x)
516      REAL bessj1,x
517      REAL ax,xx,z
518      DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
519     *                 s1,s2,s3,s4,s5,s6,y
520      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
521     *     s5,s6
522      DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
523     * 242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2,
524     * s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0,
525     * 99447.43394d0,376.9991397d0,1.d0/
526      DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
527     * .2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
528     * -.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
529      if(abs(x).lt.8.)then
530        y=x**2
531        bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
532     *         y*(s4+y*(s5+y*s6)))))
533      else
534        ax=abs(x)
535        z=8./ax
536        y=z**2
537        xx=ax-2.356194491
538        bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
539     *    p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.,x)
540      endif
541      END
542C