| 1 | C PROGRAM QUASI_LINES_2D |
|---|
| 2 | C----------------------------------------------------------------------- |
|---|
| 3 | C This is a 2-D version of QUASI_LINES, a 1-D Bayesian Quasi-elastic |
|---|
| 4 | C line-fitting program. The resolution function is given on a constant |
|---|
| 5 | C X-binning grid, and is assumed to be invariant. The background is |
|---|
| 6 | C assumed to be linear and the spectrum a sum of Laurentzians, centered |
|---|
| 7 | C at the "origin". The maximum number of lines allowed to be fit is 3. |
|---|
| 8 | C The data are assumed to be of intermediate Genie binary format, such |
|---|
| 9 | C as that output by batch-mode ICON. |
|---|
| 10 | C----------------------------------------------------------------------- |
|---|
| 11 | C Written by: D.S. Sivia, Rutherford Appleton Lab., Oxfordshire, England. |
|---|
| 12 | C----------------------------------------------------------------------- |
|---|
| 13 | C Initial release for trials by CJC, MAA, WSH. (DSS: 27-FEB-1992) |
|---|
| 14 | C Changed file access from sequential to direct. (DSS: 5-MAR-1992) |
|---|
| 15 | C Changed STEPSZ reduction from 0.85 to 0.6 . (DSS: 7-MAR-1992) |
|---|
| 16 | C Can read in a detector normalisation file. (DSS: 12-MAR-1992) |
|---|
| 17 | C Readonly for resol. file; silent-running modified. (DSS: 14-MAY-1992) |
|---|
| 18 | C WBF array-size increased in subroutine INBUFF. (DSS: 28-MAY-1992) |
|---|
| 19 | C Fixed bug in CLRBUF; problem for lots of detectors.(DSS: 29-MAY-1992) |
|---|
| 20 | C Hard-wired meV as input and ueV as output. (DSS: 2-OCT-1992) |
|---|
| 21 | C Put in STEPSZ safety-valve in case Chisq stuck! (DSS: 27-OCT-1992) |
|---|
| 22 | C Took out small and negative number safety-valve. (DSS: 17-DEC-1992) |
|---|
| 23 | C Fixed bug in DATIN1 (no more large error-bars!). (DSS: 27-SEP-1993) |
|---|
| 24 | C----------------------------------------------------------------------- |
|---|
| 25 | C Option : elastic peak , o_el=0 (no), =1 (yes) |
|---|
| 26 | C : background , o_bgd=2 (sloping), =1 (flat), =0 (zero) |
|---|
| 27 | C : width1 , o_w1=1 (fix), =0 (free) |
|---|
| 28 | C----------------------------------------------------------------------- |
|---|
| 29 | C |
|---|
| 30 | SUBROUTINE QLwat(numb,x_in,y_in,e_in,reals,opft, |
|---|
| 31 | 1 XD_in,XB_in,YB_in,Wy_in,We_in,dtn,xsc,sfile,rfile,l_fn, |
|---|
| 32 | 2 nd_out,xout,yout,eout,yfit,yprob) |
|---|
| 33 | INCLUDE 'mod_files.f90' |
|---|
| 34 | INCLUDE 'mod_data.f90' |
|---|
| 35 | INCLUDE 'options.f90' |
|---|
| 36 | COMMON /FFTCOM/ FRES(m_d2),FWRK(m_d2),XJ(m_d),TWOPIK(m_d1),NFFT |
|---|
| 37 | COMMON /DATCOM/ XDAT(m_d),DAT(m_d),SIG(m_d),NDAT |
|---|
| 38 | COMMON /DINTRP/ IPDAT(m_d),XPDAT(m_d) |
|---|
| 39 | COMMON /FITCOM/ FIT(m_d),RESID(m_d),NFEW,FITP(m_p),EXPF(m_d1,6) |
|---|
| 40 | COMMON /SCLCOM/ BSCL,ASCL,WSCL,SCLVEC(m_p,2),GSCL |
|---|
| 41 | COMMON /GRDCOM/ DDDPAR(m_d,m_p),FR2PIK(m_d2,2) |
|---|
| 42 | COMMON /QW1COM/ QW1(m_sp),SIGQW1(m_sp),ISPEC |
|---|
| 43 | COMMON /WRKCOM/ WORK(m_d2,2) |
|---|
| 44 | COMMON /BESSJ/ AJ0,AJ1,AJ2 |
|---|
| 45 | COMMON /ModPars/ NBIN,IMIN,IMAX,RSCL,BNORM |
|---|
| 46 | real x_in(m_d), y_in(m_d), e_in(m_d) |
|---|
| 47 | cf2py intent(in) :: x_in, y_in, e_in !sample data |
|---|
| 48 | real reals(4) |
|---|
| 49 | cf2py intent(in) :: reals !real parameters |
|---|
| 50 | real XD_in(m_d), XB_in(m_d), YB_in(m_d) |
|---|
| 51 | cf2py intent(in) :: XD_in, XB_in, YB_in !sample xrange, res data (blur) |
|---|
| 52 | real Wy_in(m_sp), We_in(m_sp), dtn(m_sp), xsc(m_sp) |
|---|
| 53 | cf2py intent(in) :: Wy_in, We_in, dtn, xsc !fixed width data & res scaling |
|---|
| 54 | integer numb(9) |
|---|
| 55 | cf2py intent(in) :: numb !integer parameters |
|---|
| 56 | integer opft(4) |
|---|
| 57 | cf2py intent(in) :: opft !options parameters |
|---|
| 58 | integer l_fn |
|---|
| 59 | cf2py intent(in) :: l_fn !length of filenames |
|---|
| 60 | character*140 sfile, rfile |
|---|
| 61 | cf2py intent(in) :: sfile, rfile !sample & res filenames |
|---|
| 62 | integer nd_out |
|---|
| 63 | cf2py intent(out) :: nd_out !number of output points |
|---|
| 64 | real xout(m_d), yout(m_d), eout(m_d) |
|---|
| 65 | cf2py intent(out) :: xout, yout, eout !data values |
|---|
| 66 | real yfit(4*m_d) |
|---|
| 67 | cf2py intent(out) :: yfit !fit values |
|---|
| 68 | real yprob(4) |
|---|
| 69 | cf2py intent(out) :: yprob !probability values |
|---|
| 70 | integer l_title,l_user |
|---|
| 71 | character*80 user,title |
|---|
| 72 | real XBLR(m_d),YBLR(m_d) |
|---|
| 73 | REAL GRAD(m_p),HESS(m_p,m_p),DPAR(m_p),COVAR(m_p,m_p) |
|---|
| 74 | REAL SIGPAR(m_p),FITPSV(m_p),DTNORM(m_sp),XSCALE(m_sp) |
|---|
| 75 | REAL PRBSV(4,m_sp),POUT(4,m_sp), |
|---|
| 76 | 1 PRMSV(7,4,m_sp),SIGSV(7,4,m_sp) |
|---|
| 77 | REAL AJ0Q(m_sp),AJ1Q(m_sp),AJ2Q(m_sp) |
|---|
| 78 | INTEGER INDX(m_p) |
|---|
| 79 | DATA AROT /0.98/ |
|---|
| 80 | prog='w' |
|---|
| 81 | CALL init_paras |
|---|
| 82 | do n=1,m_p |
|---|
| 83 | GRAD(n)=0.0 |
|---|
| 84 | DPAR(n)=0.0 |
|---|
| 85 | SIGPAR(n)=0.0 |
|---|
| 86 | FITPSV(n)=0.0 |
|---|
| 87 | do m=1,m_p |
|---|
| 88 | HESS(n,m)=0.0 |
|---|
| 89 | COVAR(n,m)=0.0 |
|---|
| 90 | end do |
|---|
| 91 | end do |
|---|
| 92 | do n=1,m_sp |
|---|
| 93 | do n1=1,4 |
|---|
| 94 | PRBSV(n1,n)=0.0 |
|---|
| 95 | POUT(n1,n)=0.0 |
|---|
| 96 | do n2=1,7 |
|---|
| 97 | PRMSV(n2,n1,n)=0.0 |
|---|
| 98 | SIGSV(n2,n1,n)=0.0 |
|---|
| 99 | end do |
|---|
| 100 | end do |
|---|
| 101 | end do |
|---|
| 102 | c numb = [ngrp, nsp, ntc, Ndat, nbin, Imin, Imax, Nb, nrbin] |
|---|
| 103 | NSPEC=numb(1) !no. of groups |
|---|
| 104 | ISP=numb(2) !group number |
|---|
| 105 | ISPEC=ISP |
|---|
| 106 | ntc=numb(3) !no. of points |
|---|
| 107 | NDAT=numb(4) |
|---|
| 108 | NBIN=numb(5) |
|---|
| 109 | IMIN=numb(6) |
|---|
| 110 | IMAX=numb(7) |
|---|
| 111 | NB=numb(8) |
|---|
| 112 | nrbin=numb(9) |
|---|
| 113 | c reals = [efix, theta[isp], rscl, bnorm] |
|---|
| 114 | efix=reals(1) |
|---|
| 115 | theta(ISP)=reals(2) |
|---|
| 116 | RSCL=reals(3) |
|---|
| 117 | BNORM=reals(4) |
|---|
| 118 | do n=1,m_d |
|---|
| 119 | xin(n)=x_in(n) |
|---|
| 120 | yin(n)=y_in(n) |
|---|
| 121 | ein(n)=e_in(n) |
|---|
| 122 | XDAT(n)=XD_in(n) |
|---|
| 123 | XBLR(n)=XB_in(n) |
|---|
| 124 | YBLR(n)=YB_in(n) |
|---|
| 125 | end do |
|---|
| 126 | o_el=opft(1) |
|---|
| 127 | o_bgd=opft(2) |
|---|
| 128 | o_w1=opft(3) |
|---|
| 129 | do n=1,m_sp |
|---|
| 130 | DTNORM(n)=dtn(n) !DTNORM |
|---|
| 131 | XSCALE(n)=xsc(n) !XSCALE |
|---|
| 132 | yprob=0.0 |
|---|
| 133 | end do |
|---|
| 134 | lptfile='' |
|---|
| 135 | fileout1='' |
|---|
| 136 | l_lpt=l_fn+8 |
|---|
| 137 | lptfile(1:l_lpt)=sfile(1:l_fn)//'_QLw.lpt' |
|---|
| 138 | l_file=l_fn+8 |
|---|
| 139 | fileout1(1:l_lpt)=sfile(1:l_fn)//'_QLw.ql1' |
|---|
| 140 | l_title=9 |
|---|
| 141 | title='<unknown>' |
|---|
| 142 | l_user=9 |
|---|
| 143 | user='<unknown>' |
|---|
| 144 | if(o_w1.eq.1)then |
|---|
| 145 | do I=1,NSPEC |
|---|
| 146 | QW1(I)=Wy_in(I) |
|---|
| 147 | QW1(I)=0.5*(ABS(QW1(I))+0.00001) |
|---|
| 148 | SIGQW1(I)=We_in(I) |
|---|
| 149 | SIGQW1(I)=0.5*(ABS(SIGQW1(I))+0.00001) |
|---|
| 150 | end do |
|---|
| 151 | endif |
|---|
| 152 | if(ISP.eq.1)then !print info |
|---|
| 153 | call open_f(53,lptfile) |
|---|
| 154 | write(53,1107)sfile |
|---|
| 155 | 1107 format(' Sample file : ',a140) |
|---|
| 156 | write(53,1108)rfile |
|---|
| 157 | 1108 format(' Resolution file : ',a140) |
|---|
| 158 | write(53,1110)xin(imin),xin(imax) |
|---|
| 159 | 1110 format(' Energy range: ',f8.3,' to ',f8.3,' meV') |
|---|
| 160 | if (o_el.eq.0) write(53,1111) |
|---|
| 161 | if (o_el.eq.1) write(53,1112) |
|---|
| 162 | if (o_bgd.eq.2) write(53,1113) |
|---|
| 163 | if (o_bgd.eq.1) write(53,1114) |
|---|
| 164 | if (o_bgd.eq.0) write(53,1115) |
|---|
| 165 | 1111 format(' Elastic option : NO peak') |
|---|
| 166 | 1112 format(' Elastic option : WITH peak') |
|---|
| 167 | 1113 format(' Background option : sloping') |
|---|
| 168 | 1114 format(' Background option : flat') |
|---|
| 169 | 1115 format(' Background option : zero') |
|---|
| 170 | if(o_w1.eq.1)then |
|---|
| 171 | write(53,1116) |
|---|
| 172 | else |
|---|
| 173 | write(53,1117) |
|---|
| 174 | endif |
|---|
| 175 | 1116 format(' Width option : fixed from file ') |
|---|
| 176 | 1117 format(' Width option : free') |
|---|
| 177 | close(unit=53) |
|---|
| 178 | call open_f(1,fileout1) |
|---|
| 179 | do n=1,1 |
|---|
| 180 | write(n,1107)sfile |
|---|
| 181 | write(n,1121)title(1:l_title) |
|---|
| 182 | 1121 format(' Title : ',a) |
|---|
| 183 | write(n,1122)user(1:l_user) |
|---|
| 184 | 1122 format(' User : ',a) |
|---|
| 185 | write(n,1123)NSPEC,NDAT,xin(imin),xin(imax) |
|---|
| 186 | 1123 FORMAT(2X,2I10,2x,2f10.3) |
|---|
| 187 | write(n,*)' -------------------------------------------------' |
|---|
| 188 | write(n,1108)rfile |
|---|
| 189 | write(n,1124)n |
|---|
| 190 | 1124 FORMAT(i3) |
|---|
| 191 | write(n,*)' -------------------------------------------------' |
|---|
| 192 | close(unit=n) |
|---|
| 193 | end do |
|---|
| 194 | endif !end print |
|---|
| 195 | CALL BLRINT(XBLR,YBLR,NB,1.0,.FALSE.) |
|---|
| 196 | NDAT=ntc-1 |
|---|
| 197 | CALL DPINIT |
|---|
| 198 | CALL GDINIT |
|---|
| 199 | NDAT=ntc-1 |
|---|
| 200 | CALL DATIN(ISP,DTNORM,IDUF) |
|---|
| 201 | if(o_w1.eq.1)then |
|---|
| 202 | OPEN(UNIT=53,FILE=lptfile,STATUS='old',FORM='formatted', |
|---|
| 203 | 1 access='append') |
|---|
| 204 | write(53,1120)QW1(isp) |
|---|
| 205 | 1120 format(' qlm> width 1 fixed at ',f10.5) |
|---|
| 206 | close(unit=53) |
|---|
| 207 | endif |
|---|
| 208 | CALL BLRINT(XBLR,YBLR,NB,XSCALE(ISP),.TRUE.) |
|---|
| 209 | CALL DPINIT |
|---|
| 210 | CALL PRINIT(FITP,3,NFEW,1) |
|---|
| 211 | CALL FileInit(1,ISP) |
|---|
| 212 | CALL BESINT(AJ0Q(ISP),AJ1Q(ISP),AJ2Q(ISP),AROT*QAVRG(ISP)) |
|---|
| 213 | AJ0=AJ0Q(ISP) |
|---|
| 214 | AJ1=AJ1Q(ISP) |
|---|
| 215 | AJ2=AJ2Q(ISP) |
|---|
| 216 | CALL REFINAw(GRAD,HESS,DPAR,3+NFEW,DETLOG,INDX,COVAR) |
|---|
| 217 | goto 2 |
|---|
| 218 | 1 CALL SEARCHw(GRAD,HESS,DPAR,NFEW,INDX,COVAR,FITP) |
|---|
| 219 | 2 NPARMS=4+3*NFEW |
|---|
| 220 | CHIOLD=CCHI(FITP) |
|---|
| 221 | CALL VCOPY(FITP,FITPSV,NPARMS) |
|---|
| 222 | STEPSZ=0.3 |
|---|
| 223 | IF (NFEW.GT.1) STEPSZ=STEPSZ/10.0 |
|---|
| 224 | IAGAIN=0 |
|---|
| 225 | CDIFMN=0.003 |
|---|
| 226 | do I=1,200 |
|---|
| 227 | CALL REFINEw(GRAD,HESS,NPARMS,DETLOG,INDX,COVAR,STEPSZ) |
|---|
| 228 | CALL NEWEST(COVAR,GRAD,NPARMS,NFEW,DPAR,FITP) |
|---|
| 229 | CNORM=CCHI(FITP) |
|---|
| 230 | IF (CNORM.LE.CHIOLD) THEN |
|---|
| 231 | CHIDIF=(CHIOLD-CNORM)/CNORM |
|---|
| 232 | IF (ABS(CHIDIF).LE.CDIFMN) THEN |
|---|
| 233 | IF (IAGAIN.EQ.0) THEN |
|---|
| 234 | CDIFMN=0.00005 |
|---|
| 235 | STEPSZ=0.15 |
|---|
| 236 | IAGAIN=1 |
|---|
| 237 | ELSE |
|---|
| 238 | GOTO 3 |
|---|
| 239 | ENDIF |
|---|
| 240 | ENDIF |
|---|
| 241 | CHIOLD=CNORM |
|---|
| 242 | CALL VCOPY(FITP,FITPSV,NPARMS) |
|---|
| 243 | ELSE |
|---|
| 244 | CALL VCOPY(FITPSV,FITP,NPARMS) |
|---|
| 245 | STEPSZ=STEPSZ*0.6 |
|---|
| 246 | IF (STEPSZ.LT.1.0E-10) GOTO 3 |
|---|
| 247 | ENDIF |
|---|
| 248 | end do |
|---|
| 249 | 3 continue |
|---|
| 250 | CALL REFINEw(GRAD,HESS,NPARMS,DETLOG,INDX,COVAR,0.7) |
|---|
| 251 | CALL ERRBARw(COVAR,NP,SIGPAR,COV12) |
|---|
| 252 | CALL SEEFITw(SIGPAR,CNORM,PRMSV(1,NFEW+1,ISP),SIGSV(1,NFEW+1,ISP)) |
|---|
| 253 | CALL OUTPRM(FITP,COVAR,NPARMS,NFEW,CNORM) |
|---|
| 254 | CALL REFINEw(GRAD,HESS,NPARMS,DETLOG,INDX,COVAR,0.25) |
|---|
| 255 | CALL PROBNw(CNORM,NDAT,DETLOG,NFEW,PROBLG) |
|---|
| 256 | noff=NDAT*NFEW |
|---|
| 257 | do n=1,NDAT |
|---|
| 258 | yfit(noff+n)=FIT(n) |
|---|
| 259 | end do |
|---|
| 260 | NFEW=NFEW+1 |
|---|
| 261 | if (o_el.eq.0) then !no elastic peak |
|---|
| 262 | FITP(3)=0.0 |
|---|
| 263 | FITPSV(3)=0.0 |
|---|
| 264 | endif |
|---|
| 265 | IF (NFEW.LE.1) GOTO 1 |
|---|
| 266 | nd_out=NDAT |
|---|
| 267 | do n=1,nd_out |
|---|
| 268 | xout(n)=XDAT(n) |
|---|
| 269 | yout(n)=DAT(n) |
|---|
| 270 | if(SIG(n).lt.1.0e-10)then |
|---|
| 271 | eout(n)=0.0 |
|---|
| 272 | else |
|---|
| 273 | eout(n)=SQRT(2.0/SIG(n)) |
|---|
| 274 | endif |
|---|
| 275 | end do |
|---|
| 276 | CALL PRBOUT(PRBSV,1,ISP,POUT) |
|---|
| 277 | do n=1,2 |
|---|
| 278 | yprob(n)=POUT(n,isp) |
|---|
| 279 | end do |
|---|
| 280 | kill=0 |
|---|
| 281 | RETURN |
|---|
| 282 | END |
|---|
