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 |
---|
49 | C ------------------------------------ |
---|
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 |
---|
73 | c |
---|
74 | C***<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) |
---|
131 | C 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 |
---|
146 | C |
---|
147 | C**<search for one more & refine amplitudes>**************************** |
---|
148 | C |
---|
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 |
---|
182 | C**<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 |
---|
213 | C ------------------------------------------------------- |
---|
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(*) |
---|
228 | C |
---|
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 |
---|
310 | C***<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 |
---|
360 | C ---------------------------------- |
---|
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 |
---|
369 | C ---------------------------------------------------------- |
---|
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 |
---|
424 | C ------------------------------------------------------------------ |
---|
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 |
---|
444 | C --------------------------------- |
---|
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 |
---|
462 | C --------------------------------- |
---|
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 |
---|
473 | C --------------------------------- |
---|
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 |
---|
486 | C***<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 |
---|
514 | C ------------------ |
---|
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 |
---|
542 | C |
---|