c **** x4fti used in FEVAL - lx1 set to LT/2 **** c KLUDGE in "ne funct1" to set x7 x8 ==0 c cxxx has kludge for x5 > 0 - not used c care ifail=2 in refits - change use 17-02-03 ** version that uses fitex.h to set parameter sizes c**************************************************************** c Start of General correlated fit code here c**************************************************************** C FIT T CORRELATED DATA WITH PARAMETERS X c written Chris Michael 6-7-1994 (uses NAG - library 18) C C INPUT: DAT(it,ip,nl) storage size itst npst nlst c Fits itmin to itmax (Convention 0,..,itmax) for c operators ip=1,..npth. There are nlat data samples (assumed c statistically independent) c c Correlated chi**2 fits to DAT using Inverse Correlation c CINV as evaluated by CORFIT. The Fit function (to DAT) is c obtained from FEVAL - given NPAR parameters X c Initial values for X are needed as input c c FIT details - see FEVAL - c This version uses c nexp exponentials, npth operators c if npar=nexp*(npth+1) then m(nexp) free (idv=1 in TFIT) c else if npar=nexp*(npth)-1 then m(nexp) fixed to xdv input c ( if idv=0 in TFIT) c LT needed in FEVAL (for cosh fit - ifexp gives cosh/exp choice) c see /factf/ for factorising information in fit c c iprin (print out amount) (none if iprin=0) c c CORRELATION CHOICES: icorrt, icorrp (see CORFIT prologue) c idelt ievals xlmin c c OUTPUT: Final X() parameters, printout, chi2 c c CALLS c CORFIT deals with correlations in data c (calls CDISPLAY, INVERC, EIGEN which call F01ABF and F02ABF) c (also calls if icorrt=2 COR5, FMIN2, E04FDF, LSFUN1) c c NAG minimisation routines (E04JAF) c (calls FUNCT1, FEVAL) c Fit function is defined in FEVAL c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FITEX(dat,itmin,itmax,npth,x,npars,chi2) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (npst=npsth,itst=itsth) c ** fit workspace - see dispx, feval also for x size PARAMETER ( LIW=NPARSv+4, LW=NPARSv*(NPARSv+13) ) double precision W(LW),xl(nparsv),xu(nparsv),xx(nparsv) &,xorig(nparsv) INTEGER IW(LIW) c***input from calling program***** DIMENSION dat(0:itst,npst,nlst),x(npars) COMMON/PARMS/ xdv,xlmin,npar,nexp,npthv,lt, & icorrt,icorrp,idelt,ievals,itminv,itmaxv,nlat,iprin common/factf/nfac1,nfac2,ip1v(npst),ip2v(npst),ifs(npst) c***local storage COMMON/correl/FD(0:itst,npst),CINV(0:itst,0:itst,npst,npst) DIMENSION cdum(0:itst,0:itst,npst,npst) & ,c(0:itst,0:itst,npst,npst),fita(0:itst,npst) c Data input verify if(nparsv.lt.np*(nexp+1) ) stop ' nparsv too small' if(icorrt.eq.2.and.idelt.le.1) stop 'idelt bad choice' if(npth.gt.npst) stop 'fitex storage error: npst' if(itmin.lt.0.or.itmax.lt.0) stop 'fitex error: it < 0' if(itmin.gt.itst.or.itmax.gt.itst) & stop 'fitex storage error: itst' if(itmin.gt.itmax) stop 'fitex error: itmax' ITminv=ITmin ITmaxv=ITmax npthv=npth DO 111 It1=itmin,itmax do 111 ip1=1,npth FD(It1,ip1)=0.0 DO 111 It2=itmin,itmax do 111 ip2=1,npth CDUM(It1,It2,ip1,ip2)=0.0 111 CONTINUE XNL=1.0/FLOAT(NLAT) do 100 nl=1,nlat DO 70 Ip1=1,npth do 71 it1=itmin,itmax FD(It1,ip1)=FD(It1,ip1)+DAT(it1,ip1,nl)*XNL DO 72 ip2=1,npth do 73 it2=itmin,itmax CDUM(It1,It2,ip1,ip2)=CDUM(It1,It2,ip1,ip2) & +dat(it1,ip1,nl)*dat(it2,ip2,nl)*XNL 73 continue 72 CONTINUE 71 continue 70 CONTINUE 100 CONTINUE c assume data set independent in NL (no autocorrelation) xf=1.0/float(nlat-1) DO 80 Ip1=1,npth do 81 it1=itmin,itmax DO 82 ip2=1,npth do 83 it2=itmin,itmax C(It1,It2,ip1,ip2)=(CDUM(It1,It2,ip1,ip2)- & fd(it1,ip1)*fd(it2,ip2))*xf 83 continue 82 CONTINUE 81 continue 80 CONTINUE C c is data correl matrix call corfit(c,itmin,itmax,itst, & npth,npst,cinv,icorrt,icorrp,idelt,ievals,xlmin,iprin) C C MINIMISATION FAIRLY CONTINUOUS DERIVATIVES C C CARE X(N) GETS OVERWRITTEN IF NPARS>N>NPAR ! do ix=1,npar xorig(ix)=x(ix) enddo c ifail=0 is OK but ifail=1 suppresses error massages IFAIL = 1 C JAF: LIW.GT.N+2 LW.GT.MAX(13,12*N+N(N-1)/2) (KCF less) C IBOUND = 0 NEEDS MAX AND MIN FOR ALL X IBOUND = 1 cXXX c Ibound=0 c Do ipdx=1,npar c xl(ipdx)=-1000.0 c xu(ipdx)=1000.0 c Enddo c xl(5)=0.0 c xu(5)=2.0 ndx=0 CALL E04JAF( NPAR-ndx,IBOUND, XL,XU,X(1+ndx), & FSUMSQ, IW,LIW,W,LW, IFAIL) c ifail=1 wrong inputs c 2 400*n iterations and not converged c 3 no minimum found - lowest value given c 4 overflow c 5-8 doubt if true minimum (5 likely,...8 unlikely) IF (IFAIL.GT.0) THEN WRITE(6,992) IFAIL,fsumsq,(x(i),i=1,npar) 992 format(' first attempt, ifail=',i2,' chi2=',f15.5, & /,15f6.3/,15f6.3/,15f6.3) if(iprin.gt.0) call dispx(x,xx,npars) ifail=1 CALL E04JAF( NPAR-ndx,IBOUND, XL,XU,X(1+ndx), & FSUMSQ, IW,LIW,W,LW, IFAIL) c 17-2-03 swap here to keep orig if ifail=2 was .eq. 2 if(ifail.eq.4)then do ix=1,npar x(ix)=xorig(ix) enddo endif if(iprin.gt.0) write(6,933)ifail, fsumsq 933 format(' redone fit with new ifail=',i2,' chi2=',f15.5, & ' using orig x') ENDIF chi2=fsumsq if(iprin.gt.0)then C PRINTOUT HERE (care in calling prgm too) call dispx(x,xx,npars) ndata=(itmax-itmin+1)*npth WRITE(6,99) FSUMSQ,ITmin,itmax,npar,ndata,nlat 99 FORMAT(' FULL CHI2=',F13.3,' IT1=',I2,' it2=',i2,' npars=',i2, &' ndata=',i3,' nconfigs=',i4/ &' IP1 IP2 IT FIT DATA ERROR CHI ') call feval(fita,itst,Npst,x,npar) DO 31 IP=1,npth DO 32 IT=ITmin,ITmax fit=fita(it,ip) ERR=SQRT(C(It,it,ip,ip)) CHI=(FIT-FD(It,ip))/ERR WRITE(6,98) IP1v(ip),IP2v(ip),IT,FIT,FD(It,ip),ERR,CHI 32 CONTINUE 31 CONTINUE 98 FORMAT(2I3,I8,4F13.6) endif RETURN END C****************************************************** c c written Chris Michael 16-3-94 to 6-7-94 c for general description see C. Michael and A. McKerrell July 1994 c c input corr(it1,it2,ip1,ip2) it=itmin,..itmax ip=1,..npth c c output cinv(it1,it2,ip1,ip2) according to models: c NOTE not all elements set to zero (or symmetry cases filled) c c flags icorrt icorrp c 0 0 diagonal inversion c 1 0 NN model inverse in t, diag in p c 1 2 NN model inverse in t, correl in p c 0 2 general correl in p, diag in t c 2 0 5-diag inverse in t, diag in p c 2 2 NN model inverse in t, correl in p c 3 0 inversion (ievals eigenmodes) in t, diag in p c 3 2 inversion (ievals eigenmodes) in t, correl in p c c iprin level of printout 0-8 (use 0 for bootstrap calls) c (eigen has iprin 2: evals, 6: evects, 4: display trunc M ) c (cdisplay has iprin>4 : off diagonal averages) c c idelt abs(it1-it2) used for fitting fall-off of corrt=1 c for icorrt=2 (5-diag) fit 1 and idelt c ievals number of exact eigenvalues for approximate inverse c remainder of eigenvalues modified, see below: c Use ievals=full size for exact inverse c If ievals=99 fix minimum eigenvalue to xlmin c If ievals<0 use SVD inverse with |ievals| exact terms c If ievals>0 use ievals exact eigenmodes and remainder c same average eigenvalue: returned as xlmin c xlmin minimum eigenvalue of Correlation matrix (used if c ievals =99) c c CALLS: cdisplay inverc eigen c which in turn call NAG routines to INVERT and c find eigenvalues of a real symmetric matrix c c (also calls if icorrt=2 COR5, FMIN2, E04FDF, LSFUN1) c c c SENSIBLE CHOICES c c icorrt=0 icorrp=0 is uncorrelated fit c c 1 0 fits 1 exponential to correlation in t c for each ip. Use idelt >0 and ~(itmax-itmin)/4 c c 2 0 fits 2 exponential to corelation in t c for each ip. Use idelt >1 and ~(itmax-itmin)/4 c c 3 0 uses ievals eigenmodes of correlation in t c exactly for each ip. Use 2< ievals <= itmax-itmin+1 c Larger ievals are OK if more data points - c typically use ievals < sqrt (nlat) c Full exact inversion of correlation in t: c ievals=itmax-itmin+1 c c 3 2 uses ievals eigenmodes of correlation between c data at all t and ip. Use ievals >>2 c <= np*(itmax-itmin+1) c Larger ievals are OK if more data points - c typically use ievals < sqrt (nlat) c Full exact inversion of correlation in t and p: c ievals=np*(itmax-itmin+1) c C****************************************************** subroutine corfit(corr,itmin,itmax,itst, & npth,npst,cinv,icorrt,icorrp,idelt,ievals,xlmin,iprin) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' parameter (npthl=npsth,itstl=itsth*npthl) dimension corr(0:itst,0:itst,npst,npst), & cinv(0:itst,0:itst,npst,npst), & cnorm(0:itstl),cdum(0:itstl,0:itstl), & xas(npthl),xrr(npthl,npthl),xl(0:itstl) c Data input verify if(itstl.lt.itmax.or.npthl.lt.npth) & stop 'corfit storage error' itval=itmax-itmin itvalp=itval+1 itstp=itstl+1 if(iprin.ge.4)then C DISPLAY input correlation blocks do 50 ip1=1,npth do 51 ip2=ip1,npth do 56 it1=itmin,itmax do 57 it2=itmin,itmax cdum(it1,it2)=corr(it1,it2,ip1,ip2) & /sqrt( corr(it1,it1,ip1,ip1)* corr(it2,it2,ip2,ip2) ) 57 continue 56 continue write(6,910)ip1,ip2 910 format(/' input normalised correlation t*t for paths=',2i2) call cdisplay(cdum,itmin,itmax,itstl,5) c evaluate estimate of smallest eigenvalue from alternating eigenmode v=0.0 s1=-1.0 do 65 it1=itmin,itmax s1=-s1 s2=-1.0 do 66 it2=itmin,itmax s2=-s2 v=v+s1*cdum(it1,it2)*s2 66 continue 65 continue c normalise vector +-+-+-,.. v=v/itvalp write(6,971)v 971 format(/' Estimates of smallest eigenvalues by various methods' & /' projection of norm correl on vector +-+-,.. is',f14.7) c evaluate estimates of smallest eigenvector from submatrices if(ip1.eq.ip2.and.icorrt.ge.2)then c maximum submatrix given by statistics: ievals? idtmax=min(itval,ievals) do 68 itd=1,idtmax itdiff=itd+1 itb=itmax-itd do 67 ita=itmin,itb c use cinv as dummy to receive output not needed c use iprin=-1 specially to return least eigenvalue c (no output to cinv) call eigen(cdum(ita,ita),itstp,itdiff,1,xl(ita),cinv,-1) 67 continue write(6,939)itdiff,itdiff 939 format(i2,' by ',i2,' submatrix has min evals:') write(6,9)(xl(it),it=itmin,itb) 9 format(5f13.6) 68 continue endif 51 continue 50 continue endif if(icorrp.eq.0.or.npth.eq.1)then c**************BLOCK DIAGONAL IN PATHS **************************** do 1 ip=1,npth if(icorrt.eq.0) then c---diagonal in t ------------------------------------ do 20 it=itmin,itmax cinv(it,it,ip,ip)=1.0/corr(it,it,ip,ip) 20 continue else c measure corr normalised by diagonal: model it c icorrt.ne.0 ------------------------------------ do 30 it=itmin,itmax cnorm(it)=1.0/sqrt(corr(it,it,ip,ip)) 30 continue calculate average normalised correlation at idelt if(idelt+itmin.gt.itmax)stop ' idelt too big' nrav=0 rav=0.0 nrav1=0 rav1=0.0 do 301 it1=itmin,itmax it2=it1+idelt if(it2.le.itmax)then nrav=nrav+1 rav=rav+corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2) endif it2=it1+1 if(it2.le.itmax)then nrav1=nrav1+1 rav1=rav1+corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2) endif 301 continue rav=rav/float(nrav) rav1=rav1/float(nrav1) if(iprin.ge.2)then write(6,930)ip,1,rav1,idelt,rav 930 format(' norm correl ip=',i2,' average C(t1-t2=',i2,')=',f8.5 & /' average C(t1-t2=',i2,')=',f8.5) endif if(icorrt.eq.1)then c---nearest neighbour in cinv in t----------- c... option to guess NN model to try and fit C ravv=(rav)**(1.0/float(idelt)) xdiag0=1.0/(1.0-ravv*ravv) xdiag=(1.0+ravv*ravv)*xdiag0 xnn=-ravv*xdiag0 c speed up by just diagonal and off diagonal initialised (see FUNCT1) do 430 it1=itmin,itmax-1 it2=it1+1 cinv(it1,it1,ip,ip)=xdiag*(cnorm(it1)*cnorm(it1)) cinv(it1,it2,ip,ip)=xnn*(cnorm(it1)*cnorm(it2)) cinv(it2,it1,ip,ip)=xnn*(cnorm(it1)*cnorm(it2)) 430 continue cinv(itmin,itmin,ip,ip)=xdiag0*(cnorm(itmin)*cnorm(itmin)) cinv(itmax,itmax,ip,ip)=xdiag0*(cnorm(itmax)*cnorm(itmax)) if(iprin.ge.4) then do 410 it1=itmin,itmax do 411 it2=itmin,itmax cdum(it1,it2)=0.0 if(abs(it1-it2).eq.1) cdum(it1,it2)=xnn if(it1.eq.it2) cdum(it1,it2)=xdiag if(it1.eq.itmax.and.it2.eq.itmax)cdum(it1,it2)=xdiag0 if(it1.eq.itmin.and.it2.eq.itmin)cdum(it1,it2)=xdiag0 cinv(it1,it2,ip,ip)=cdum(it1,it2)*(cnorm(it1)*cnorm(it2)) 411 continue 410 continue if(iprin.ge.8) then call cdisplay(cdum,itmin,itmax,itstl,iprin) c to get more info on eigenvectors: call eigen(cdum(itmin,itmin),itstp,itvalp, & itvalp,xlmin,cdum(itmin,itmin),2) else call inverc(cdum(itmin,itmin),cdum(itmin,itmin), & itstp,itvalp) endif write(6,911)ip 911 format(/' model normalised correlation t*t for path=',i2) call cdisplay(cdum,itmin,itmax,itstl,5) endif endif if(icorrt.eq.2)then call cor5(rav1,1,rav,idelt,dd0,xn0,dd1,dd,xn,xnn,iprin) c speed up by just diagonal and 4 off diagonals initialised (see FUNCT1) c could speed up more by just upper diagonals used (symmetric) if(iprin.ge.4)then do 510 it1=itmin,itmax do 511 it2=itmin,itmax cinv(it1,it2,ip,ip)=0.0 511 continue 510 continue endif if(itmax-itmin.lt.4) stop ' 5-diagonal wont fit' cinv(itmin,itmin,ip,ip)=dd0*(cnorm(itmin)*cnorm(itmin)) cinv(itmin+1,itmin+1,ip,ip)=dd1*(cnorm(itmin+1)*cnorm(itmin+1)) do 500 it1=itmin+2,itmax-2 cinv(it1,it1,ip,ip)=dd*(cnorm(it1)*cnorm(it1)) 500 continue cinv(itmax-1,itmax-1,ip,ip)=dd1*(cnorm(itmax-1)*cnorm(itmax-1)) cinv(itmax,itmax,ip,ip)=dd0*(cnorm(itmax)*cnorm(itmax)) cinv(itmin,itmin+1,ip,ip)=xn0*(cnorm(itmin)*cnorm(itmin+1)) cinv(itmin+1,itmin,ip,ip)=xn0*(cnorm(itmin)*cnorm(itmin+1)) do 530 it1=itmin+1,itmax-2 it2=it1+1 cinv(it1,it2,ip,ip)=xn*(cnorm(it1)*cnorm(it2)) cinv(it2,it1,ip,ip)=xn*(cnorm(it1)*cnorm(it2)) 530 continue cinv(itmax,itmax-1,ip,ip)=xn0*(cnorm(itmax)*cnorm(itmax-1)) cinv(itmax-1,itmax,ip,ip)=xn0*(cnorm(itmax)*cnorm(itmax-1)) do 501 it1=itmin,itmax-2 it2=it1+2 cinv(it1,it2,ip,ip)=xnn*(cnorm(it1)*cnorm(it2)) cinv(it2,it1,ip,ip)=xnn*(cnorm(it1)*cnorm(it2)) 501 continue if(iprin.ge.4)then do 520 it1=itmin,itmax do 521 it2=itmin,itmax cdum(it1,it2)=cinv(it1,it2,ip,ip)/(cnorm(it1)*cnorm(it2)) 521 continue 520 continue if(iprin.ge.8) then write(6,913)ip 913 format(/' model norm inverse correlation t*t for path=',i2) call cdisplay(cdum,itmin,itmax,itstl,iprin) c to get more info on eigenvectors: call eigen(cdum(itmin,itmin),itstp,itvalp, & itvalp,xlmin,cdum(itmin,itmin),2) else call inverc(cdum(itmin,itmin),cdum(itmin,itmin), & itstp,itvalp) endif write(6,912)ip 912 format(/' model normalised correlation t*t for path=',i2) call cdisplay(cdum,itmin,itmax,itstl,5) endif endif c-----expand Corr in eigenvalues, and modify if needed--------------- if(icorrt.eq.3)then do 610 it1=itmin,itmax do 611 it2=itmin,itmax cdum(it1,it2)=corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2) 611 continue 610 continue c inverse (truncated if ievals relevant) call eigen(cdum(itmin,itmin),itstp,itvalp, & ievals,xlmin,cdum(itmin,itmin),min(iprin,2)) do 620 it1=itmin,itmax do 621 it2=itmin,itmax cinv(it1,it2,ip,ip)=cdum(it1,it2)*(cnorm(it1)*cnorm(it2)) 621 continue 620 continue if(iprin.ge.8)then write(6,915)ip 915 format(/' inverse normalised correlation t*t for path=',i2) call cdisplay(cdum,itmin,itmax,itstl,5) endif if(iprin.ge.4) then write(6,919)ip 919 format(/' model normalised correlation t*t for path=',i2) call inverc(cdum(itmin,itmin),cdum(itmin,itmin),itstp,itvalp) call cdisplay(cdum,itmin,itmax,itstl,5) endif endif endif 1 continue RETURN else if(icorrp.ge.2.and.npth.gt.1)then c**************OFF DIAGONAL SYMMETRIC IN PATHS ******************** if(icorrt.eq.0)then if(iprin.ge.8) write(6,977) 977 format(/' it inverse norm corr(ip1,ip2) ,....') do 700 it=itmin,itmax do 702 ip1=1,npth do 704 ip2=1,npth cdum(ip1,ip2)=corr(it,it,ip1,ip2) 704 continue 702 continue call inverc(cdum(1,1),cdum(1,1),itstp,npth) do 712 ip1=1,npth do 714 ip2=1,npth cinv(it,it,ip1,ip2)=cdum(ip1,ip2) if(iprin.ge.8) cdum(ip1,ip2)=cdum(ip1,ip2)* & sqrt(corr(it,it,ip1,ip1)*corr(it,it,ip2,ip2)) 714 continue 712 continue if(iprin.ge.8)then write(6,978)it, ((cdum(i,j),j=1,npth),i=1,npth) 978 format(i4,3(6f12.5/)) endif 700 continue return endif c hereon ip and it correlated: if(itstl.lt.npth*itvalp) stop 'corfit: storage 2,2 problem' do 732 ip1=1,npth do 730 it1=itmin,itmax iv1=(ip1-1)*itvalp+it1-itmin cnorm(iv1)=1.0/sqrt(corr(it1,it1,ip1,ip1)) 730 continue 732 continue if(icorrt.eq.1.or.icorrt.eq.2)then c model diagonal and off diagonal correlations calculate average normalised correlation at idelt (and at 0 ) c positivity of model needed: one way is to fix slope as shown do 735 ip1=1,npth do 734 ip2=ip1,npth nrav=0 rav=0.0 rrr=0.0 do 736 it1=itmin,itmax iv1=(ip1-1)*itvalp+it1-itmin iw2=(ip2-1)*itvalp+it1-itmin if(ip1.ne.ip2) & rrr=rrr+corr(it1,it1,ip1,ip2)*cnorm(iv1)*cnorm(iw2) it2=it1+idelt iv2=(ip2-1)*itvalp+it2-itmin if(it2.le.itmax)then nrav=nrav+1 rav=rav+corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2) endif it2=it1-idelt iv2=(ip2-1)*itvalp+it2-itmin if(it2.ge.itmin)then nrav=nrav+1 rav=rav+corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2) endif 736 continue rav=rav/float(nrav) if(ip1.ne.ip2) then xrr(ip1,ip2)=rrr/itvalp else c****icorrt=1 ****here***** xas(ip1)=-log(rav)/float(idelt) xrr(ip1,ip2)=1.0 endif rav=rav/xrr(ip1,ip2) if(iprin.ge.2)then write(6,931)ip1,ip2,idelt,rav,xrr(ip1,ip2) 931 format(' norm correl ip1,ip2=',2i2,' C(',i2, & ')/C(0)=',f8.5,' C(0)=',f8.5) endif 734 continue 735 continue c model correlation here do 748 ip1=1,npth do 749 ip2=ip1,npth xa=0.5*(xas(ip1)+xas(ip2)) do 740 it1=itmin,itmax iv1=(ip1-1)*itvalp+it1-itmin ivx1=(ip2-1)*itvalp+it1-itmin do 741 it2=itmin,itmax iv2=(ip2-1)*itvalp+it2-itmin ivx2=(ip1-1)*itvalp+it2-itmin cdum(iv1,iv2)=xrr(ip1,ip2)*exp(-xa*abs(it2-it1)) c assume symmetric in ip1, ip2 cdum(ivx1,ivx2)=cdum(iv1,iv2) 741 continue 740 continue 749 continue 748 continue c cdum set up icort = 1 or 2 endif if (icorrt.eq.3)then c**************FULL MATRIX IN TIME AND PATHS ************************** do 800 it1=itmin,itmax do 801 it2=itmin,itmax do 802 ip1=1,npth iv1=(ip1-1)*itvalp+it1-itmin do 804 ip2=1,npth iv2=(ip2-1)*itvalp+it2-itmin cdum(iv1,iv2)=corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2) 804 continue 802 continue 801 continue 800 continue endif c invert cdum for icorrt = 1,2 or 3 c need full inversion even for icorrt=1 because of off diag in ip care negative eigenvalues possible if model is dud c call inverc(cdum,cdum,itstp,npth*itvalp) iev=ievals if(icorrt.eq.1) iev=npth*itvalp call eigen(cdum,itstp,npth*itvalp,iev, & xlmin,cdum,min(2,iprin)) do 810 it1=itmin,itmax do 811 it2=itmin,itmax do 812 ip1=1,npth iv1=(ip1-1)*itvalp+it1-itmin do 814 ip2=1,npth iv2=(ip2-1)*itvalp+it2-itmin cinv(it1,it2,ip1,ip2)=cdum(iv1,iv2)*cnorm(iv1)*cnorm(iv2) 814 continue 812 continue 811 continue 810 continue c display inverse model correlation (normalised) if(iprin.ge.8)then do 870 ip1=1,npth do 871 ip2=ip1,npth write(6,985)ip1,ip2 985 format(/' inverse norm correl for ip=',i2,',',i2) call cdisplay(cdum((ip1-1)*itvalp,(ip2-1)*itvalp) & ,0,itval,itstl,iprin) 871 continue 870 continue endif c display truncated input - reinvert cdum (not needed icorrt=1) if(iprin.ge.4.and.iev.ne.npth*itvalp)then call inverc(cdum,cdum,itstp,npth*itvalp) do 873 ip1=1,npth do 874 ip2=ip1,npth write(6,988)ip1,ip2 988 format(/' model norm correl for ip=',i2,',',i2) call cdisplay(cdum((ip1-1)*itvalp,(ip2-1)*itvalp) & ,0,itval,itstl,iprin) 874 continue 873 continue endif RETURN else stop ' corfit: wrong icorrt icorrp choices' endif return end c****************************************************************** c CALLED by CORFIT c Code to find sensible 5-diagonal matrix inverse c to given (reduced) correlations c C(idt=it1)=p1, C(idt=it2)=p2 (assumes C(0)=1.0 ) c c CALLS FMIN2 c c 5-diagonal convention: c [x y c 0 0 0 0 c y z b c 0 0 0 c c b a b c 0 0 c 0 c b a b c 0 c 0 0 c b a b c c 0 0 0 c b z y c 0 0 0 0 c y x] c c****************************************************************** subroutine cor5(p1in,it1in,p2in,it2in,x,y,z,a,b,c,iprin) implicit double precision (a-h,o-z) common/abfit/p1,p2,it1,it2 if(p1.gt.1.0.or.p2.gt.1.0) stop 'cor5: corr > 1' p1=p1in p2=p2in it1=it1in it2=it2in c assume 1 exponential to start xl1=(p2)**(1.0/float(it2)) call fmin2(xl1,0.0,xf,yf,iprin) xl1=xf xl2=yf if(abs(xl1-xl2).lt.1.0D-10) xl2=xl1-1.0D-6 do 10 ir=1,19 ctw=( xl1**(ir+1)*(1.0-xl2**2) -xl2**(ir+1)*(1.0-xl1**2) ) & /( (xl1-xl2)*(1.0+xl1*xl2) ) if(ir.eq.it1)c1v=ctw if(ir.eq.it2)c2v=ctw if(ir.eq.1)p1v=ctw if(ir.eq.2)p2v=ctw if(iprin.ge.8) write(6,900)ir,ctw 900 format(' idt=',i2,' corr~ =',f12.6) 10 continue xd=1.0/( (xl1-xl2)*(1.0+xl1*xl2) ) xc1=xd*(1.0-xl2**2) xc2=xd*(1.0-xl1**2) FFF=(xc1*xl1**(it1+1)-xc2*xl2**(it1+1)-p1)**2 FFF=FFF+(xc1*xl1**(it2+1)-xc2*xl2**(it2+1)-p2)**2 if(iprin.ge.4) write(6,92)xl1,xl2,FFF,it1,c1v,it2,c2v 92 format(' Eigenvalues=',2g14.6,' chi^2 fit=',g15.7 / & ' with c(',i2,')=',f12.6,' and c(',i2,')=',f12.6) if(abs(min(xl1,xl2)).le.1.0D-10) then xl=max(xl1,xl2) xx1=xl+1.0/xl bbya=-1.0/xx1 a=1.0/(1.0+p1v*2.0*bbya) b=bbya*a c=0.0 x=a/(xx1*xl) z=a y=b else xx1=xl1+1.0/xl1 xx2=xl2+1.0/xl2 bbyc=-xx1-xx2 cbya=1.0/(2.0+xx1*xx2) bbya=bbyc*cbya a=1.0/(1.0+p1v*2.0*bbya+p2v*2.0*cbya) b=bbya*a c=cbya*a x=c/(xl1*xl2) z=a-c*xl1*xl2 y=-c*(xl1+xl2)/(xl1*xl2) endif if(iprin.ge.8) write(6,1)p1v,p2v,x,y,z,a,b,c 1 format(2f12.5,3x,3f15.5/3f15.5) return end c**************************************************************** c CALLED by COR5 c NAG minimising code F(x,y) start x0,y0 final xf,yf c LEAST SQUARES ( tested 30-6-94 ) c c THIS VERSION FIXES 0.0 <= x,y <= 1.0 c by mapping xl= tanh(x)*.5+.5 c CALLS LSFUN1, E04FDF c**************************************************************** subroutine fmin2(x0,y0,xf,yf,iprin) implicit double precision (a-h,o-z) external LSFUN1 dimension xv(2),W(50),iw(50) c xv(1)=x0 c xv(2)=y0 xv(1)=0.5 * log ( x0 /(1.0-x0) ) xv(2)=-1.0 ifail=1 call e04fdf(2,2,xv,fsumsq,iw,50,w,50,ifail) if(ifail.gt.1) write(6,99) ifail, fsumsq 99 format('cor5: icorrt=2 fit ifail=',i8,' chi2=',g15.6) xf=tanh(xv(1))*0.5+0.5 yf=tanh(xv(2))*0.5+0.5 return end c****************************************************************** c called by e04fdf c i=1,m Sum_of_f(i)**2 fitted parameters xm(n) c****************************************************************** subroutine lsfun1(m,n,xm,f) implicit double precision (a-h,o-z) dimension xm(n),f(m),x(4) common/abfit/p1,p2,it1,it2 x(1)=tanh(xm(1))*0.5+0.5 x(2)=tanh(xm(2))*0.5+0.5 xd=1.0/( (x(1)-x(2))*(1.0+x(1)*x(2)) ) xc1=xd*(1.0-x(2)**2) xc2=xd*(1.0-x(1)**2) F(1)=xc1*x(1)**(it1+1)-xc2*x(2)**(it1+1)-p1 F(2)=xc1*x(1)**(it2+1)-xc2*x(2)**(it2+1)-p2 return end c****************************************************************** c CALLED by CORFIT c display c(it1,it2) it=itmin,..itmax and calculate average c off axis diagonal values (if iprin > 4) c c normalised for display by scale XSC c c c****************************************************************** subroutine cdisplay(c,itmin,itmax,itst,iprin) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' parameter (itstv=npsth*itsth) dimension c(0:itst,0:itst),irv(itstv),inum(0:itstv),rav(0:itstv) if(itmax.gt.itstv)stop 'cdisplay: storage error' nit=itmax-itmin do 1 it=0,itmax inum(it)=it 1 continue c fix SCALE to display as I3 format xsc=100.0 if(abs(c(itmin,itmin)).gt.1.01) xsc=10.0 if(abs(c(itmin,itmin)).gt.10.01) xsc=1.0 if(abs(c(itmin,itmin)).gt.100.1) xsc=0.1 if(abs(c(itmin,itmin)).gt.1000.1) xsc=0.01 write(6,991)itmin,itmax,itmin,itmax,xsc 991 format(' cdisplay: c(',i2,':',i2,',',i2,':',i2,') with scale' & ,' xsc=',f8.3) write(6,993) itst,(inum(it),it=itmin,itmax) 993 format(i4,':',25i3) write(6,994) 994 format(' ') do 2 idt=0,nit rav(idt)=0.0 2 continue do 31 it1=itmin,itmax do 32 it2=itmin,itmax idt=abs(it2-it1) rav(idt)=rav(idt)+c(it1,it2) irv(it2)=nint(xsc*c(it1,it2)) 32 continue write(6,99) it1,(irv(it2),it2=itmin,itmax) 99 format(i4,':',25i3) 31 continue if(iprin.gt.4)then write(6,907) 907 format(' Average off-diagonal elements') do 3 idt=1,nit rav(idt)=0.5*rav(idt)/(nit-idt+1) 3 continue rav(0)=rav(0)/(nit+1) write(6,902)(inum(it),rav(it),it=0,nit) 902 format(5(i3,':',f11.5)) endif return end C======================================================================C C C SUBROUTINE FUNCT1( NPAR , X, FC ) C C C LEAST SQUARES FUNCTION CALLED BY E04JAF C Calls Feval to get F(X); C needs data FD and correlation CINV c needs flags icorrt icorrp C CALLS FEVAL C======================================================================C C IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' parameter (itst=itsth,npst=npsth) c***input from calling program***** DIMENSION X(NPAR) COMMON/correl/FD(0:itst,npst),CINV(0:itst,0:itst,npst,npst) COMMON/PARMS/xdv,xlmin,nparv,nexp,npth,lt, & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin c**local*** Logical ifp Dimension FIT(0:itst,npst) ifp=icorrp.ne.0.and.npth.gt.1 call Feval(fit,itst,npst,X,npar) FC = 0. c g4 gi gi g5 gi - want a1 zero to g4 gi - say 7,8 c KLUDGE here to force a parameter to fixed value (assumed zero here) if(ikludge .eq. 123) FC=100000.0*( x(7)**2+x(8)**2 ) if(ikludge .eq. 246) FC=100000.0*( x(5)**2+x(6)**2 ) if(ikludge .eq. 987) FC=100000.0*( x(9)**2+x(10)**2 ) npthm=npth-1 c itminl=9 itminll=10 from fitex.h do 21 ip1=1,npth DO 20 it1 = itmin,itmax ifdo=1 if(ikludge .eq. 666 .and. ip1 .gt. 10 .and. it1 .lt. itminl) & ifdo=0 if(ikludge .eq. 666 .and. (ip1 .gt. 19 .or. ip1 .eq. 15) & .and. it1 .lt. itminll) ifdo=0 if (ifdo .eq. 1) & fc=fc+cinv(it1,it1,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))**2 20 continue 21 continue if(.not.ifp)then if(icorrt.eq.0) return c---Nearest neighbour in t (assume symmetric)------ if(icorrt.ge.1) then do 31 ip1=1,npth DO 30 it1 = itmin,itmax-1 it2=it1+1 fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))* & 2.0*(fit(it2,ip1)-fd(it2,ip1)) 30 continue 31 continue endif if(icorrt.ge.2)then c also add further diagonal in t (assume symmetric) do 300 ip1=1,npth DO 303 it1 = itmin,itmax-2 it2=it1+2 fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))* & 2.0* (fit(it2,ip1)-fd(it2,ip1)) 303 continue 300 continue endif c---General correlation in t (assume symmetric)------ if(icorrt.gt.2.and.(itmax-itmin).gt.2) then do 41 ip1=1,npth DO 40 it1 = itmin,itmax-3 DO 42 it2 = it1+3,itmax fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))* & 2.0* (fit(it2,ip1)-fd(it2,ip1)) 42 continue 40 continue 41 continue endif else c+++++++several paths, correlated here ++++++ c add to diagonal correlation term first if(icorrt.eq.0)then c offdiagonal in p do 50 ip1=1,npthm do 51 ip2=ip1+1,npth do 55 it=itmin,itmax fc=fc+cinv(it,it,ip1,ip2)*(fit(it,ip1)-fd(it,ip1))* & 2.0*(fit(it,ip2)-fd(it,ip2)) 55 continue 51 continue 50 continue else c offdiagonal in t do 471 ip1=1,npth DO 470 it1 = itmin,itmax-1 DO 472 it2 = it1+1,itmax fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))* & 2.0* (fit(it2,ip1)-fd(it2,ip1)) 472 continue 470 continue 471 continue c offdiagonal in p c assume symmetric under ip1 ip2 swap not it1 it2 do 400 ip1=1,npthm DO 401 ip2 =ip1+1,npth DO 402 it1 = itmin,itmax do 403 it2 = itmin,itmax fc=fc+cinv(it1,it2,ip1,ip2)*(fit(it1,ip1)-fd(it1,ip1))* & 2.0* (fit(it2,ip2)-fd(it2,ip2)) 403 continue 402 continue 401 continue 400 continue endif endif if(FC.lt.-1.0D-8) STOP ' Funct1: chi**2 is negative' RETURN END c****************************************************************** c CALLED by FUNCT1, FITEX, TFIT, GRAPH c Function for fits: F(it,ip) it convention 0,.. c given NPAR parameters X c as written NEXP exponential fit to NPTH operators ip c options: cosh (needs lt ) ifexp=0 c exp ifexp=-1 (Could speed up this case in code) c sinh (needs lt) ifexp=1 c m(nexp) fixed at xdv if npar < nexp*(npth+1) c c Conventions if factorising: c F=x1p*x2p exp(-xm1*t)+ x1p*x2p exp(-(xm1+xm2)*t) c else (nfac1=1) c F=x2p exp(-xm1*t)+ x2p exp(-(xm1+xm2)*t) c c ORDER x path1 m0; path2 m0 , etc c normalise coefficients to 1.0 at t=tmin for ip=1 (by XD in TFIT) c c changes here need initial values set appropriately (TFIT FITEX) c and display modified (DISPX) c npar set by funct1 if called E04 c****************************************************************** subroutine feval(f,itst,npst,xv,npar) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (npstv=npsth) dimension xm(nparsv),f(0:itst,npst),xv(npar),x(nparsv) c***input from calling program***** COMMON/PARMS/ xdv,xlmin,nparv,nexp,npth,lt, & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin common/factf/nfac1,nfac2,ip1v(npstv),ip2v(npstv),ifs(npstv) np=max(nfac1,nfac2) nm=nexp*np nmm=nm+nexp do ix=1,npar x(ix)=xv(ix) enddo xm(1)=x(nm+1) if( ifs(1).eq. 88) then do ip = 1,np do it=itmin,itmax f(it,ip)=xm(1)+x(ip)*(float(it-lt/2))**2 enddo enddo return endif do 10 it=itmin,itmax ca=exp(-xm(1)*(it-itmin)) cb= exp(-xm(1)*(lt-it-itmin)) cc=ca+cb cs=ca-cb do 11 ip=1,npth if(ifs(ip).eq.0)then ce=cc elseif(ifs(ip).eq.1)then ce=cs elseif(ifs(ip).eq.-1)then ce=ca elseif(ifs(ip).eq.9)then c point-point (code in tmfit24 - so care if not linked) irv=2 c lx1=24 lx1=lt/2 c use interpolation here now (need lx1=lt/2) -- code in tmfit48 only c BUT STUB HERE ALSO FOR COMPILER ce=x4fti(xm(1),it,irv,lx1,lt)/x4fti(xm(1),itmin,irv,lx1,lt) endif if(nfac1.eq.1)then f(it,ip)=x(ip2v(ip))*ce else f(it,ip)=x(ip1v(ip))*x(ip2v(ip))*ce endif 11 continue 10 continue if(nexp.gt.1)then do 20 im=2,nexp xm(im)=x(nm+im) 20 continue if(npar.eq.nmm-1) xm(nexp)=xdv do 30 it=itmin,itmax do 33 im=2,nexp ca=exp(-xm(im)*(it-itmin)) cb= exp(-xm(im)*(lt-it-itmin)) cc=ca+cb cs=ca-cb do 31 ip=1,npth if(ifs(ip).eq.0)then ce=cc elseif(ifs(ip).eq.1)then ce=cs elseif(ifs(ip).eq.-1)then ce=ca elseif(ifs(ip).eq.9)then c point-point (code in tmfit24 - so care if not linked) irv=2 c lx1=24 lx1=lt/2 ce=x4fti(xm(im),it,irv,lx1,lt)/x4fti(xm(im),itmin,irv,lx1,lt) endif if(nfac1.eq.1)then f(it,ip)=f(it,ip)+x(np*(im-1)+ip2v(ip))*ce else f(it,ip)=f(it,ip)+x(np*(im-1)+ip1v(ip))* & x(np*(im-1)+ip2v(ip))*ce endif 31 continue 33 continue 30 continue endif Return end c-------------------------------------------------------------------- c CALLED by CORFIT c invert C (in place allowed) dimension ndstv size ndata c assumes C is positive definite and symmetric c cdum has only i>=j elements with ci(i,j) filled c CALLS F01ABF c-------------------------------------------------------------------- subroutine inverc(c,ci,ndstv,ndata) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (npst=npsth,itst=itsth,NDST=NPST*ITST) dimension c(ndstv,ndata),ci(ndstv,ndata) double precision WKSPCE(NDST),cdum(ndst,ndst) if(ndst.lt.ndata.or.ndstv.lt.ndata) stop ' inverc: storage error' IFAIL=0 CALL F01ABF(c,NDSTv,NDATA,cdum,NDST,WKSPCE,IFAIL) if(ifail.eq.1) write(6,98)ifail 98 format(' inversion in inverc failed ifail=',i2,' not positive & definite') if(ifail.eq.2) write(6,99)ifail 99 format(' inversion in inverc failed ifail=',i2,' ill conditioned') do 1 i=1,ndata do 2 j=1,i ci(i,j)=cdum(i,j) ci(j,i)=cdum(i,j) 2 continue 1 continue RETURN END C*********************************************************************** c CALLED by CORFIT C CALCULATE EIGENVECTORS and hence INVERSE c If ieval=99 fix minimum eigenvalue to xlmin c If ieval<0 use SVD inverse keping |ieval| exact terms c If ieval>0 use ieval exact eigenmodes and remainder c same average eigenvalue: returned as xlmin c c output cinv (with modified eigenmodes) can be in place c c NOTE iprin=-1 returns least eigenvalue as xlmin c iprin>0 negative eigenvalue warning and inversion summary c iprin>1 eigenvalues and truncation info c iprin>3 display (truncated) input matrix c iprin>7 eigenvectors c c MODS 2-3-94 to "smooth" xl more c CALLS cdisplay, f02abf C*********************************************************************** SUBROUTINE EIGEN(C,NPVST,NPV,ieval,xlmin,cinv,iprin) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (NPSTS=itsth*npsth) DIMENSION C(NPVST,NPVst),cinv(npvst,npvst) Dimension EVECTS(NPSTS,NPSTS),xev(NPSTS), & TEMP(NPSTS,NPSTS),EVALS(NPSTS),work(NPSTS) ievals=iabs(ieval) if(NPSTS.lt.npv.or.npvst.lt.npv) stop 'Eigen: storage error' if(ievals.gt.npv.and.ieval.ne.99) stop 'Eigen: ievals too big' DO 10 I=1,NPV DO 10 J=1,NPV 10 TEMP(I,J)=C(I,J) IFAIL=0 c Real symmetric (uses lower triangle): eigenvalues ascending order c and normalised eigenvectors Call f02abf(temp,NPSTS,npv,evals,evects,NPSTS,work,ifail) if(ifail.ne.0.and.iprin.gt.0)write(6,919)ifail 919 format(' error in Eigen: Ifail=',i2) if(iprin.eq.-1) then xlmin=evals(1) return endif if(evals(1).le.0.0.and.iprin.gt.0)then write(6,95)npv,evals(1) 95 format(' CARE: ',i2,'th eigenvalue =',f12.5,' NEGATIVE') endif if(iprin.gt.1)write(6,981) 981 format(' eigenvalues ') IF(IPRIN.GT.1)WRITE(6,98)NPV,(EVALS(IP),IP=1,NPV) if(iprin.gt.7)write(6,982) 982 format(' eigenvectors ') do 483 i=1,npv IF(IPRIN.GT.7)WRITE(6,98)I,(EVECTS(IP,I),IP=1,NPV) 483 continue 98 FORMAT(I6,6F12.6,7(/,6X,6f12.6)) do 200 ip1=1,npv do 201 ip2=1,npv temp(ip1,ip2)=0.0 cinv(ip1,ip2)=0.0 201 continue 200 continue if(ieval.lt.npv.and.ieval.ge.0)then xevs=0.0 do 150 ie=1,npv-ievals xevs=xevs+evals(ie) 150 continue xlmin=xevs/(npv-ievals) if(iprin.ge.1) write(6,97)npv-ievals,evals(1), & evals(npv-ievals),xlmin 97 format(' Inversion:',i3,' smallest eigenvalues ',f10.6, & ' to', f10.6,/' replaced by max(value, average): ', f10.7) if(xevs.le.0.0.and.iprin.ge.1)write(6,902) 902 format(' warning: average eigenvalue zero !') else if(ieval.eq.npv) then c full inversion if(iprin.ge.1) write(6,92) evals(1) 92 format(' Full inversion: smallest eigenvalue=',f10.7) xlmin=0.0 else if(ieval.eq.99)then c ievals < 0 so just use xlmin here if(iprin.ge.1) write(6,93) xlmin 93 format(' Inversion: eigenvalues <',f10.6,' replaced by ', & 'that value') else if(ieval.lt.0)then if(iprin.ge.1) write(6,94)npv-ievals,evals(1),evals(npv-ievals) 94 format(' SVD inverse: ',i2,' smallest eigenvalues ',f10.6, &' to ',f10.6,/' replaced by zero') else stop ' eigen: wrong ievals selection' endif c modify eigenvales as needed xnsum=0.0 xosum=0.0 do 155 ie=1,npv xl=evals(ie) c replace xl < xlmin by xlmin if(ieval.eq.99.and.xl.lt.xlmin) xl=xlmin c replace npv-ievals smallest evals by average if less if(ieval.ge.0.and.ie.le.(npv-ievals).and.xl.lt.xlmin) & xl=xlmin c replace npv-ievals smallest eigenvalues by SVD (zero in inverse) if(ieval.lt.0.and.ie.le.(npv-ievals)) xl=0.0 xev(ie)=xl xnsum=xnsum+xl xosum=xosum+evals(ie) 155 continue c sum of eigenvalues of xnorm=xosum/xnsum if(xosum.ne.xnsum.and.iprin.ge.1) write(6,909) & (xev(ie)*xnorm,ie=1,npv) 909 format(' New eigenvalues: ',7(/,6X,6f12.6) ) do 100 ie=1,npv do 110 ip1=1,npv do 120 ip2=1,npv xl=xev(ie)*xnorm if(xl.le.0.0) then c replace npv-ievals smallest eigenvalues by SVD (zero in inverse) xlinv=0.0 else xlinv=1.0/xl endif temp(ip1,ip2)=temp(ip1,ip2)+ & evects(ip1,ie)*evects(ip2,ie)*xl cinv(ip1,ip2)=cinv(ip1,ip2)+ & evects(ip1,ie)*evects(ip2,ie)*xlinv 120 continue 110 continue 100 continue npstm=NPSTS-1 npvm=npv-1 if(iprin.ge.4) call cdisplay(temp,0,npvm,npstm,5) RETURN END c******************************************************** c c FT (on lattice) of 1/(p^2+m^2) c c also include form-factor exp(-a p^2) - but should not be necessary c c********************************************************* function x4ft(xmpi,it,ir,l,lx4) implicit double precision (a-h,o-z) twopi=8.0d0*atan(1.0d0) anorm=1.0d0/dfloat(l*l*l*lx4) c=0.0 lx4h=lx4/2 lh=l/2 xpf=twopi/dfloat(l) do itv=0,lx4h xtct=2.0d0 if(itv .eq. 0 .or. itv .eq. lx4h) xtct=1.0d0 xpt=itv*twopi/dfloat(lx4) xden=xmpi*xmpi+(2.0d0*sin(xpt*0.5d0) )**2 c cos is even - so fold do ixv=0,lh xxct=2.0d0 if(ixv .eq. 0 .or. ixv .eq. lh) xxct=1.0d0 xpx=ixv*xpf xden1=xden+(2.0d0*sin(xpx*0.5d0) )**2 ccv=xtct*xxct*cos(xpt*dfloat(it))*cos(xpx*dfloat(ir)) c fold to speed up sin vanishes at 0,lh - so double [was -lh,lh-1 ] do iyv=0,lh xyct=2.0d0 if(iyv .eq. 0 .or. iyv .eq. lh) xyct=1.0d0 xpy=iyv*xpf xden2=xden1+(2.0d0*sin(xpy*0.5d0) )**2 do izv=0,lh xzct=2.0d0 if(izv .eq. 0 .or. izv .eq. lh) xzct=1.0d0 xpz=izv*xpf xden3=xden2+(2.0d0*sin(xpz*0.5d0) )**2 c=c+ xyct*xzct*ccv /xden3 Enddo Enddo Enddo Enddo x4ft=c*anorm return End cccccccccccccccccccccccccccccccccc c c interpolated x4ft -- for speed c cccccccccccccccccccccccccccccccccc function x4fti(xm,itv,ir,lx1,lx4) implicit double precision (a-h,o-z) integer it,im,ir,itv,init,nmass,im0,im1,itmax,lx1,lx4 parameter (nmass=200, xm0=0.05, xmax=0.55, itmax=24) dimension xval(nmass,itmax) save xval, init data init /0/ dm=(xmax-xm0)/float(nmass-1) if (init .eq. 0 ) then init=1 do im=1,nmass xmpi=xm0+(im-1)*dm do it=1,itmax val= x4ft(xmpi,it,ir,lx1,lx4) xtr=sqrt(float(it*it+ir*ir) ) c ad hoc interpolation formula xnaive=xtr**(-1.5)*exp(-xmpi*xtr)*sqrt(xmpi) valrat=val/xnaive c write(6,90)it,xmpi,val,valrat 90 format(i3,f12.5,e15.5,f12.5) xval(im,it)=valrat enddo enddo endif c find grid points and linear interpolate if inside box xtr=sqrt(float(itv*itv+ir*ir) ) xmint = (xm-xm0)/dm if( xmint .lt. nmass .and. xmint .gt. 0 .and. & itv .le. itmax .and. itv .gt. 0) then im0=xmint+1 xmp0=xm0+(im0-1)*dm xml=xval(im0,itv)*xtr**(-1.5)*exp(-xmp0*xtr)*sqrt(xmp0) im1=im0+1 xmp1=xm0+(im1-1)*dm xmu=xval(im1,itv)*xtr**(-1.5)*exp(-xmp1*xtr)*sqrt(xmp1) c write(*,*) xmp0,xm,xmp1 c write(*,*) xml,(dxm/dm),xmu dxm=(xmint-float(im0-1))*dm val=(xml*(dm-dxm)+xmu*(dxm) ) /dm c valx4= x4ft(xm,itv,ir,lx1,lx4) c write(6,91)itv,xm,xval,valx4 91 format(i3,3g15.6) else val=x4ft(xm,itv,ir,lx1,lx4) endif x4fti=val return End