program freezer c c Finds the best hydrodynamical parameters to describe one-particle c hadron distributions and two-particle correlation data for the c model described in LA-UR-96-0782, also in http://xxx.lanl.gov c with designation nucl-th/9603007. c c INPUT FILES c freezer.in gives options and starting guesses for parameters c pim.dat pi- invariant distribution data (optional) c pip.dat pi+ invariant distribution data (optional) c kp.dat K+ invariant distribution data (optional) c km.dat K- invariant distribution data (optional) c pro.dat proton invariant distribution data (optional) c pbar.dat antiproton invariant distribution data (optional) c pim.cor pi- 2-particle correlation data (optional) c pip.cor pi+ 2-particle correlation data (optional) c kp.cor K+ 2-particle correlation data (optional) c km.cor K- 2-particle correlation data (optional) c c OUTPUT TO SCREEN (and to freezer.log) c For each iteration: c Baryon number, isospin/baryon c Values of parameters c chi^2 for each data set c Total chi^2 c Upon reaching a minimum: c Total chi^2 and chi^2/(degree of freedom) c Best values of the parameters c Calculated values of related quantities (e.g. baryon density) c Curvature matrix and its inverse the Covariance matrix c 99% confidence errors on the parameters and related quantities c c OUTPUT TO FILES c Parameters and curvature matrix to file fort.101 (for debugging) c Calculated values of P and C at each point: c pimdist.calc pi- inv. dist. c pipdist.calc pi+ inv. dist. c kpdist.calc K+ inv. dist. c kmdist.calc K- inv. dist. c prodist.calc proton inv. dist. c pbardist.calc pbar inv. dist. c pimcor.calc pi- correlation c pipcor.calc pi+ correlation c kpcor.calc K+ correlation c kmcor.calc K- correlation c IMPLICIT NONE COMMON /int/ n1,n2 COMMON /asym/ Zp,Ap,Ztar,Atar,yp,yt COMMON /simplex/chi2min,param_min c INTEGER NP,Npar,Npoints,ndata(10),i,Itmax,Ndim,option(4) PARAMETER (NP=17,Npoints=1000,Itmax=100,Ndim=5) INTEGER mdist(6),mcorr(4),k,ndof,n1(6),n2(6),j INTEGER nrot,iflag,mask(NP) REAL T,vT,et0,muB,tau,R,aT,rweak,x(10,Ndim,Npoints),mass,gammq REAL sig(10,Npoints),param(NP),covar(NP,NP),alpha(NP,NP),chisq,dof REAL chisq0,ftol,pam(NP),vertex(NP+1,NP),alamda,dfs,dnormpi,dnormK REAL vchi2(NP+1),chi2min,chi2,ftola,lam,y(10,Npoints),ycm,fac REAL dqy,conf,sys(6),lamK,ycor,sys2(4),yP,yT,ypr,ytr,fs,normpi REAL dT,dvT,det0,dmuB,dtau,dR,daT,dycm,dlam,dlamK,calqmn,normK REAL Zp,Ap,Ztar,Atar,dummy(1),fact,eigen(NP),U(NP,NP),plus,minus,f REAL pi2,ycalc(10,Npoints),chisqi(10,Npoints),param_min(NP) REAL*4 cov(NP,NP),lowpr,lowpr2,data0(10,7,Npoints) CHARACTER*20 FILE(10),FILEOUT(10) PARAMETER(pi2=6.2831853, ftol=.0001, fac=197.33) DATA mask /17*0/ DATA dT,dvT,det0,dmuB,dtau,dR,daT,dycm,dlam,dlamK,dfs,dnormpi, * dnormK /13*0./ DATA alamda,ftola,dummy /-1.,.0001,0./ DATA FILE /'pim.dat','pip.dat','kp.dat','km.dat','pro.dat', * 'pbar.dat','pim.cor','pip.cor','kp.cor','km.cor'/ DATA FILEOUT /'pimdist.calc','pipdist.calc','kpdist.calc', * 'kmdist.calc','prodist.calc','pbardist.calc', * 'pimcor.calc','pipcor.calc','kpcor.calc','kmcor.calc'/ EXTERNAL comp,chi2 character*32 logname integer index c----------------------------------------------------------------------c open(unit=1,file='freezer.index',status='unknown') index = 0 read (1,*,err=10,end=10) index 10 index = index + 1 rewind(1) write (1,*) index close (1) write (logname,'(''freezer'',i3.3,''.log'')') index write (6,'('' opening log file '',a32,''.'')') logname open(unit=100, file=logname,status='unknown') open(unit=101, file='parameters.out',status='unknown') c c Read in parameters c OPEN (unit=10,file='freezer.in',status='old') read(10,*) Zp,Ap,Ztar,Atar,yp,yt read(10,*) T,vT,et0 read(10,*) muB,tau,R read(10,*) aT,ycm,fs read(10,*) lam,lamK,normpi,normK read(10,*) (mask(i),i=1,13) read(10,*) mdist read(10,*) sys read(10,*) mcorr read(10,*) sys2 read(10,*) ycor, rweak read(10,*) option read(10,*) n1 read(10,*) n2 CLOSE(10) c c Check to make sure parameters are in allowed ranges c if ((Zp.le.0.).or.(Ztar.le.0.).or.(Ap.le.0.).or.(Atar.le.0.)) then write(*,*) 'No antimatter: Zp,Ap,Ztar,Atar must be > 0' write(100,*) 'No antimatter: Zp,Ap,Ztar,Atar must be > 0' stop endif if ((T.le.0.).or.(et0.le.0.).or.(tau.le.0.).or.(R.le.0.).or. * (lam.le.0.).or.(lamK.le.0.).or.(normpi.le.0.).or. * (normK.le.0.)) then write(*,*) 'T,et0,tau,R,lmp,lmK,normpi,normK must all be > 0' write(100,*) 'T,et0,tau,R,lmp,lmK,normpi,normK must all be > 0' stop endif if (mub*t.ge.939.0) then write (* ,*)' MuB/T should be less than 939./T' write (100,*)' MuB/T should be less than 939./T' stop endif if ((vT.le.0.).or.(vT.ge.1.)) then if (vT.eq.0.) then write(*,*) 'Sorry, expanding sources only (vT > 0.)' write(*,*) 'Try using something like vT/c = 0.001' write(100,*) 'Sorry, expanding sources only (vT > 0.)' write(100,*) 'Try using something like vT/c = 0.001' endif write(*,*) 'vT/c must be in the range 0 < vT/c < 1' write(100,*) 'vT/c must be in the range 0 < vT/c < 1' stop endif if (abs(aT).ge.1.) then if (aT.eq.-1.) then write(*,*) 'Try using aT = -0.999' write(100,*) 'Try using aT = -0.999' endif write(*,*) 'aT must be in the range -1. < aT < 1.' write(100,*) 'aT must be in the range -1. < aT < 1.' stop endif c c Echo back input parameters c write(*,1001) Zp,Ap,Ztar,Atar,yp,yt write(100,1001) Zp,Ap,Ztar,Atar,yp,yt 1001 format('Zp,Ap,Ztar,Atar,yp,ytar',4f6.0,2f12.7) lowpr = ycor lowpr2 = rweak write(*,*) 'Ymeas, rweak = ',lowpr,lowpr2 write(100,*) 'Ymeas, rweak = ',lowpr,lowpr2 write(*,*) 'option',option write(100,*) 'option',option write(*,*) 'n1',n1 write(100,*) 'n1',n1 write(*,*) 'n2',n2 write(100,*) 'n2',n2 c c Convert from units of fm to MeV^{-1} tau = tau/fac R = R/fac c c Signal to average over the relative sign of qout and qlong if (option(4).lt.0) n2(6) = -n2(6) c c If there is no correlation data, we can't vary lam or lamK if ((mcorr(1).eq.0).and.(mcorr(2).eq.0).and.(mask(10).ne.0)) then write(*,*) 'No pion correlation data -- fixing lampi' write(100,*) 'No pion correlation data -- fixing lampi' mask(10) = 0 endif if ((mcorr(3).eq.0).and.(mcorr(4).eq.0).and.(mask(11).ne.0)) then write(*,*) 'No kaon correlation data -- fixing lamK' write(100,*) 'No kaon correlation data -- fixing lamK' mask(11) = 0 endif c c Convert input parameters to solving parameters param(1) = log(T/140.) param(2) = .5*log(vT/(1.-vT)) param(3) = log(et0/2.) param(4) = muB param(5) = log(tau/.02027) param(6) = log(R/.03041) param(7) = .5*log((1.+aT)/(1.-aT)) param(8) = ycm param(9) = log(fs) param(10)= log(lam) param(11) = log(lamK) param(12) = log(normpi) param(13) = log(normK) param(14) = rweak param(15) = 0. param(16) = 0. ypr = yp ytr = yt c c c Read from one-particle distribution data files c do j=1,6 if ((j.eq.1).or.(j.eq.2)) mass = 139.6 if ((j.eq.3).or.(j.eq.4)) mass = 494. if ((j.eq.5).or.(j.eq.6)) mass = 939. if (mdist(j).eq.1) then OPEN (unit=10,file=FILE(j),status='old') OPEN (unit=10+j,file=FILEOUT(j),status='unknown') do i=1,Npoints read(10,*,end=981) x(j,1,i),x(j,2,i),y(j,i),sig(j,i) data0(j,1,i) = x(j,1,i) data0(j,2,i) = x(j,2,i) data0(j,3,i) = y(j,i) c c Transform from GeV to MeV if (option(2).eq.1) then x(j,2,i) = 1000.*x(j,2,i) if (option(3).eq.0) then y(j,i) = y(j,i)/1000000. sig(j,i) = sig(j,i)/1000000. else y(j,i) = y(j,i)/1000. sig(j,i) = sig(j,i)/1000. endif endif c c Is x(j,2,i) pT or mT-m? if (option(3).eq.0) then x(j,2,i) = sqrt(x(j,2,i)*(x(j,2,i) + 2.*mass)) elseif (x(j,2,i).eq.0.) then write(*,*) 'pT=0 point is bad data unless option(3)=0' write(100,*)'pT=0 point is bad data unless option(3)=0' else y(j,i) = y(j,i)/(pi2*x(j,2,i)) sig(j,i) = sig(j,i)/(pi2*x(j,2,i)) endif c c Include systematic errors sig(j,i) = sqrt(sig(j,i)**2 + (sys(j)*y(j,i))**2) data0(j,4,i) = sig(j,i)*data0(j,3,i)/y(j,i) enddo write(*,*) * 'WARNING: Number of points in ',FILE(j),' > Npoints' write(100,*) * 'WARNING: Number of points in ',FILE(j),' > Npoints' 981 CLOSE (10) c c Number of data points in each file ndata(j) = i-1 lowpr = sys(j) write(*,*) ndata(j), ' data points read from ',FILE(j), * ' sys err = ',lowpr write(100,*) ndata(j), ' data points read from ',FILE(j), * ' sys err = ',lowpr endif enddo c c c Read from two-particle correlation data files c do j=7,10 if (mcorr(j-6).eq.1) then if ((j.eq.7).or.(j.eq.8)) mass = 139.6 if ((j.eq.9).or.(j.eq.10)) mass = 494. OPEN (unit=10,file=FILE(j),status='old') OPEN (unit=10+j,file=FILEOUT(j),status='unknown') do i=1,Npoints read(10,*,end=982) (x(j,k,i),k=1,5),y(j,i),sig(j,i) do k=1,5 data0(j,k,i) = x(j,k,i) enddo data0(j,6,i) = y(j,i) c c Include systematic errors sig(j,i) = sqrt(sig(j,i)**2 + (sys2(j-6)*y(j,i))**2) data0(j,7,i) = sig(j,i) c c Convert x(j,3,i) from qlong to dqy=y1-y2 if (abs(option(4)).ge.2) then call convrt(mass,ycor,x(j,1,i), * x(j,2,i),x(j,4,i),x(j,5,i),x(j,3,i),dqy,option(4)) x(j,3,i) = dqy endif enddo write(*,*) * 'WARNING: Number of points in ',FILE(j),' > Npoints' write(100,*) * 'WARNING: Number of points in ',FILE(j),' > Npoints' 982 CLOSE (10) c c Number of data points in each file ndata(j) = i-1 lowpr = sys2(j-6) write(*,*) ndata(j), ' data points read from ',FILE(j), * ' sys err = ',lowpr write(100,*) ndata(j), ' data points read from ',FILE(j), * ' sys err = ',lowpr endif enddo c c c Determine the number of parameters and the degrees of freedom c Npar = 0 do j=1,13 if (mask(j).eq.1) Npar = Npar + 1 enddo if ((Npar.eq.0).and.(option(1).ne.0)) then write(*,*) 'Must allow some parameters to vary to minimize' write(*,*) 'Setting option(1) = 0 to continue' write(100,*) 'Must allow some parameters to vary to minimize' write(100,*) 'Setting option(1) = 0 to continue' option(1) = 0 endif c ndof = 0 do k=1,10 ndof = ndof + ndata(k) enddo ndof = ndof - Npar dof = float(ndof) write(*,3141) ndof+Npar,Npar,dof write(100,3141) ndof+Npar,Npar,dof 3141 format (' Degrees of freedom =',i7,' - ',i2,' =',f8.0) c c c Minimize chi^2 by the Simplex method (not using derivatives) c if ((option(1).eq.1).or.(option(1).eq.3)) then write(*,*) write(100,*) write(*,*) '***** Minimization by the simplex method *****' write(100,*) '***** Minimization by the simplex method *****' c c Condense vector param into vector pam made up only of variables i = 0 do j=1,13 if (mask(j).eq.1) then i = i + 1 pam(i) = param(j) endif enddo c c Quick and dirty just to get close if (option(1).eq.3) then ftola = .001*float(Npar) c c k = maximum allowed number of iterations k = 40 else c c maximum allowed number of iterations is set to default in amoeba k = 0 endif c c Evaluate starting points of the simplex do i=1,Npar+1 pam(i) = pam(i) + .1 do j=1,Npar vertex(i,j) = pam(j) enddo vchi2(i) = chi2(pam,x,y,sig,ndata,param,mask) pam(i) = pam(i) - .1 enddo c c Simplex minimization routine calls the function chi2 c CALL amoeba(vertex,vchi2,NP+1,NP,Npar,ftola,chi2,k * ,x,y,sig,ndata,param,mask) c c Choose the best point found so far c do j=1,13 if (mask(j).eq.1) then param(j) = param_min(j) endif enddo endif c c c Minimization of chi^2 by Levenberg-Marquardt (using derivatives) c k = 0 if (option(1).gt.1) then write(*,*) write(100,*) write(*,*) '***** Minimization by the Levenberg-Marquardt', * ' method *****' write(100,*) '***** Minimization by the Levenberg-Marquardt', * ' method *****' c c Iterate Itmax times or until desired precision is achieved do i=1,Itmax c c Minimization routine uses the subroutine comp c CALL mrqmin(x,y,sig,ndata,param,mask,NP,covar,alpha,NP, * chisq,comp,alamda,ycalc,chisqi) c c Minimum chi^2 so far found write(*,1002) chisq write(100,1002) chisq 1002 format ('chi^2_min = ',f20.12) c c Write intermediate parameters to fort.101 write(101,*) 'iteration number ',i write(101,*) 'param' write(101,*) param write(101,*) c c Condition for considering point the minimum if ((i.gt.1).and.(abs(chisq-chisq0)/chisq.lt.ftol)) then j = j + 1 if (j.ge.8) goto 999 else j = 0 endif chisq0 = chisq c c Shutdown message when aT -> -1. if ((abs(param(7)).gt.4.).and.(mask(7).eq.1)) then k = k + 1 if ((k.ge.3).or.(abs(param(7)).gt.100.)) then write(*,*) write(100,*) if (param(7).lt.0.) then write(*,*) '*******************************************' write(*,*) 'aT wants to go to -1, but aT must be > -1' write(*,*) 'Try restarting the program with aT = -0.999' write(*,*) ' and do not allow aT to vary' write(*,*) '*******************************************' write(100,*)'*******************************************' write(100,*) 'aT wants to go to -1, but aT must be > -1' write(100,*)'Try restarting the program with aT = -0.999' write(100,*) ' and do not allow aT to vary' write(100,*)'*******************************************' else write(*,*) '*******************************************' write(*,*) 'aT wants to go to 1, but aT must be < 1' write(*,*) 'Try restarting the program with aT = 0.999' write(*,*) ' and do not allow aT to vary' write(*,*) '*******************************************' write(100,*)'*******************************************' write(100,*) 'aT wants to go to 1, but aT must be < 1' write(100,*)'Try restarting the program with aT = 0.999' write(100,*) ' and do not allow aT to vary' write(100,*)'*******************************************' endif stop endif else k = 0 endif enddo 999 continue c c Calculate curvature matrix at the minimum alamda = 0. CALL mrqmin(x,y,sig,ndata,param,mask,NP,covar,alpha,NP, * chisq,comp,alamda,ycalc,chisqi) write(101,*) 'Final parameters and curvature matrix' write(101,*) 'param' write(101,*) param write(101,*) write(101,*) 'alpha' do i=1,Npar write(101,*) (alpha(i,j),j=1,Npar) enddo endif c c Calculate the curvature matrix when not using Lev-Marq minimization c if ((option(1).eq.0).or.(option(1).eq.1)) then if (option(1).eq.0) then write(*,*) write(100,*) write(*,*) '***** Summary Run (no minimization) *****' write(100,*) '***** Summary Run (no minimization) *****' endif alamda = -2. CALL mrqmin(x,y,sig,ndata,param,mask,NP,covar,alpha,NP, * chisq,comp,alamda,ycalc,chisqi) write(*,1002) chisq write(100,1002) chisq write(101,*) 'Final parameters and curvature matrix' write(101,*) 'param' write(101,*) param write(101,*) write(101,*) 'alpha' do i=1,Npar write(101,*) (alpha(i,j),j=1,Npar) enddo endif c c Write out the calculated values of one-part dists. at each point c do j=1,6 do i=1,ndata(j) if (option(2).eq.1) then lowpr = 1000000.*ycalc(j,i) else lowpr = ycalc(j,i) endif if (option(3).ne.0) lowpr = lowpr*pi2*data0(j,2,i) lowpr2 = chisqi(j,i) write(10+j,*) (data0(j,k,i),k=1,4),lowpr,lowpr2 enddo CLOSE(10+j) enddo c c Write out the calculated values of correlations at each point c do j=7,10 do i=1,ndata(j) lowpr = ycalc(j,i) lowpr2 = chisqi(j,i) write(10+j,*) (data0(j,k,i),k=1,7),lowpr,lowpr2 enddo CLOSE(10+j) enddo c c Best values of the parameters c T = 140.*exp(param(1)) vT = .5*(1. + tanh(param(2))) et0 = 2.*exp(param(3)) muB = param(4) tau = fac*.02027*exp(param(5)) R = fac*.03041*exp(param(6)) aT = tanh(param(7)) ycm = param(8) fs = exp(param(9)) lam = exp(param(10)) lamK = exp(param(11)) normpi = exp(param(12)) normK = exp(param(13)) c c c Putting in correct jacobians for the curvature matrix alpha c j = 0 if (mask(1).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/T alpha(i,j) = alpha(i,j)/T enddo endif if (mask(2).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/(2.*vT*(1.-vT)) alpha(i,j) = alpha(i,j)/(2.*vT*(1.-vT)) enddo endif if (mask(3).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/et0 alpha(i,j) = alpha(i,j)/et0 enddo endif if (mask(4).eq.1) j = j + 1 if (mask(5).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/tau alpha(i,j) = alpha(i,j)/tau enddo endif if (mask(6).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/R alpha(i,j) = alpha(i,j)/R enddo endif if (mask(7).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/(1.-aT*aT) alpha(i,j) = alpha(i,j)/(1.-aT*aT) enddo endif if (mask(8).eq.1) j = j + 1 if (mask(9).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/fs alpha(i,j) = alpha(i,j)/fs enddo endif if (mask(10).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/lam alpha(i,j) = alpha(i,j)/lam enddo endif if (mask(11).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/lamK alpha(i,j) = alpha(i,j)/lamK enddo endif if (mask(12).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/normpi alpha(i,j) = alpha(i,j)/normpi enddo endif if (mask(13).eq.1) then j = j + 1 do i=1,Npar alpha(j,i) = alpha(j,i)/normK alpha(i,j) = alpha(i,j)/normK enddo endif c c c Invert the curvature matrix to get the covariance matrix c do i=1,Npar do j=1,Npar covar(i,j) = alpha(i,j) enddo enddo CALL gaussj(covar,Npar,NP,dummy,0,1) c c To get 90 percent confidence levels, multiply delta by sqrt(conf) c The value of conf depends on the number of parameters. c c if (Npar.eq.0) conf = 0. c if (Npar.eq.1) conf = 2.70554 c if (Npar.eq.2) conf = 4.60517 c if (Npar.eq.3) conf = 6.25139 c if (Npar.eq.4) conf = 7.77944 c if (Npar.eq.5) conf = 9.23635 c if (Npar.eq.6) conf = 10.6446 c if (Npar.eq.7) conf = 12.0170 c if (Npar.eq.8) conf = 13.3616 c if (Npar.eq.9) conf = 14.6837 c if (Npar.eq.10) conf = 15.9871 c if (Npar.eq.11) conf = 17.2750 c if (Npar.eq.12) conf = 18.5494 c if (Npar.eq.13) conf = 19.8119 c c 99 percent confidence c if (Npar.eq.0) conf = 0. if (Npar.eq.1) conf = 6.63490 if (Npar.eq.2) conf = 9.21034 if (Npar.eq.3) conf = 11.3449 if (Npar.eq.4) conf = 13.2767 if (Npar.eq.5) conf = 15.0863 if (Npar.eq.6) conf = 16.8119 if (Npar.eq.7) conf = 18.4753 if (Npar.eq.8) conf = 20.0902 if (Npar.eq.9) conf = 21.6660 if (Npar.eq.10) conf = 23.2093 if (Npar.eq.11) conf = 24.7250 if (Npar.eq.12) conf = 26.2170 if (Npar.eq.13) conf = 27.6883 c conf = sqrt(conf) c c c Diagonal elements of the covariance matrix give errors c i = 0 if (mask(1).eq.1) then i = i + 1 dT = conf*sqrt(covar(i,i)) endif if (mask(2).eq.1) then i = i + 1 dvT = conf*sqrt(covar(i,i)) endif if (mask(3).eq.1) then i = i + 1 det0 = conf*sqrt(covar(i,i)) endif if (mask(4).eq.1) then i = i + 1 dmuB = conf*sqrt(covar(i,i)) endif if (mask(5).eq.1) then i = i + 1 dtau = conf*sqrt(covar(i,i)) endif if (mask(6).eq.1) then i = i + 1 dR = conf*sqrt(covar(i,i)) endif if (mask(7).eq.1) then i = i + 1 daT = conf*sqrt(covar(i,i)) endif if (mask(8).eq.1) then i = i + 1 dycm = conf*sqrt(covar(i,i)) endif if (mask(9).eq.1) then i = i + 1 dfs = conf*sqrt(covar(i,i)) endif if (mask(10).eq.1) then i = i + 1 dlam = conf*sqrt(covar(i,i)) endif if (mask(11).eq.1) then i = i + 1 dlamK = conf*sqrt(covar(i,i)) endif if (mask(12).eq.1) then i = i + 1 dnormpi = conf*sqrt(covar(i,i)) endif if (mask(13).eq.1) then i = i + 1 dnormK = conf*sqrt(covar(i,i)) endif c write(*,*) write(100,*) write(*,*) 'Here is the BEST point !' write(100,*) 'Here is the BEST point !' write(*,*) '**********************************************' write(100,*) '**********************************************' write(*,*) write(100,*) write(*,*) 'Total chi^2 = ',chisq write(100,*) 'Total chi^2 = ',chisq write(*,*) 'chi^2/d.o.f = ',chisq/dof write(100,*) 'chi^2/d.o.f = ',chisq/dof write(*,*) 'Probability of perfect model = ', * gammq(.5*dof,.5*chisq) write(100,*) 'Probability of perfect model = ', * gammq(.5*dof,.5*chisq) c if (Npar.ne.13) then write(*,*) write(100,*) write(*,*) 'Fixed parameters:' write(100,*) 'Fixed parameters:' endif if (mask(1).eq.0) write(*,2001) T if (mask(2).eq.0) write(*,2002) vT if (mask(3).eq.0) write(*,2003) et0 if (mask(4).eq.0) write(*,2004) muB if (mask(5).eq.0) write(*,2005) tau if (mask(6).eq.0) write(*,2006) R if (mask(7).eq.0) write(*,2007) aT if (mask(8).eq.0) write(*,2008) ycm if (mask(9).eq.0) write(*,2009) fs if (mask(10).eq.0) write(*,2010) lam if (mask(11).eq.0) write(*,2011) lamK if (mask(12).eq.0) write(*,2012) normpi if (mask(13).eq.0) write(*,2013) normK if (mask(1).eq.0) write(100,2001) T if (mask(2).eq.0) write(100,2002) vT if (mask(3).eq.0) write(100,2003) et0 if (mask(4).eq.0) write(100,2004) muB if (mask(5).eq.0) write(100,2005) tau if (mask(6).eq.0) write(100,2006) R if (mask(7).eq.0) write(100,2007) aT if (mask(8).eq.0) write(100,2008) ycm if (mask(9).eq.0) write(100,2009) fs if (mask(10).eq.0) write(100,2010) lam if (mask(11).eq.0) write(100,2011) lamK if (mask(12).eq.0) write(100,2012) normpi if (mask(13).eq.0) write(100,2013) normK 2001 format (' T =', f12.7,' MeV') 2002 format (' vT/c =', f12.7) 2003 format (' et0 =', f12.7) 2004 format ('muB/T =', f12.7) 2005 format (' tau =', f12.7,' fm/c') 2006 format (' R =', f12.7,' fm') 2007 format (' aT =', f12.7) 2008 format (' ycm =', f12.7) 2009 format (' fs =', f12.7) 2010 format (' lam =', f12.7) 2011 format (' lamK =', f12.7) 2012 format ('normpi=', f12.7) 2013 format ('normK =', f12.7) c if (Npar.gt.0) then write(*,*) write(100,*) write(*,*) 'Variable parameters:' write(100,*) 'Variable parameters:' endif if (mask(1).eq.1) write(*,3001) T, dT if (mask(2).eq.1) write(*,3002) vT, dvT if (mask(3).eq.1) write(*,3003) et0, det0 if (mask(4).eq.1) write(*,3004) muB, dmuB if (mask(5).eq.1) write(*,3005) tau, dtau if (mask(6).eq.1) write(*,3006) R, dR if (mask(7).eq.1) write(*,3007) aT, daT if (mask(8).eq.1) write(*,3008) ycm, dycm if (mask(9).eq.1) write(*,3009) fs, dfs if (mask(10).eq.1) write(*,3010) lam, dlam if (mask(11).eq.1) write(*,3011) lamK, dlamK if (mask(12).eq.1) write(*,3012) normpi, dnormpi if (mask(13).eq.1) write(*,3013) normK, dnormK if (mask(1).eq.1) write(100,3001) T, dT if (mask(2).eq.1) write(100,3002) vT, dvT if (mask(3).eq.1) write(100,3003) et0, det0 if (mask(4).eq.1) write(100,3004) muB, dmuB if (mask(5).eq.1) write(100,3005) tau, dtau if (mask(6).eq.1) write(100,3006) R, dR if (mask(7).eq.1) write(100,3007) aT, daT if (mask(8).eq.1) write(100,3008) ycm, dycm if (mask(9).eq.1) write(100,3009) fs, dfs if (mask(10).eq.1) write(100,3010) lam, dlam if (mask(11).eq.1) write(100,3011) lamK, dlamK if (mask(12).eq.1) write(100,3012) normpi, dnormpi if (mask(13).eq.1) write(100,3013) normK, dnormK 3001 format (' T = ', f12.7, ' +/- ',f11.7,' MeV') 3002 format (' vT/c = ', f12.7, ' +/- ',f11.7) 3003 format (' et0 = ', f12.7, ' +/- ',f11.7) 3004 format ('muB/T = ', f12.7, ' +/- ',f11.7) 3005 format (' tau = ', f12.7, ' +/- ',f11.7,' fm/c') 3006 format (' R = ', f12.7, ' +/- ',f11.7,' fm') 3007 format (' aT = ', f12.7, ' +/- ',f11.7) 3008 format (' ycm = ', f12.7, ' +/- ',f11.7) 3009 format (' fs = ', f12.7, ' +/- ',f11.7) 3010 format (' lam = ', f12.7, ' +/- ',f11.7) 3011 format (' lamK = ', f12.7, ' +/- ',f11.7) 3012 format ('normpi= ', f12.7, ' +/- ',f11.7) 3013 format ('normK = ', f12.7, ' +/- ',f11.7) c write(*,*) write(100,*) write(*,*) 'Calculated Quantities:' write(100,*) 'Calculated Quantities:' fact = tanh(ycm) write(*,3999) fact,tanh(ycm+dycm)-fact,tanh(ycm-dycm)-fact write(100,3999) fact,tanh(ycm+dycm)-fact,tanh(ycm-dycm)-fact fact = tanh(et0) write(*,4000) fact,tanh(et0+det0)-fact,tanh(et0-det0)-fact write(100,4000) fact,tanh(et0+det0)-fact,tanh(et0-det0)-fact c c Find Eigenvalues and Eigenvectors of the matrix alpha c do i=1,Npar do j=1,Npar cov(i,j) = covar(i,j) covar(i,j) = alpha(i,j) enddo enddo if (Npar.ge.1) then CALL Jacobi(covar,Npar,NP,eigen,U,nrot) do i=1,Npar eigen(i) = sqrt(eigen(i))/conf enddo endif c c Find quantities at minimum (f) c and their extrema on the 99% confidence hyperellipse (plus,minus) c do i=1,13 iflag = i c c If aT = 0, t2-t1 must equal 0 also c if ((i.eq.4).and.(aT.eq.0.)) then f = 0. plus = 0. minus = 0. else CALL hyperr(T,vT,et0,muB,tau,R,aT,ycm,lam,lamK,U,eigen, * Npar,mask,iflag,f,plus,minus,ypr,ytr,fs) endif c if (Npar.le.0) then plus = 0. minus = 0. elseif (Npar.eq.1) then plus = calqmn(T+dT,vT+dvT,et0+det0,muB+dmuB,tau+dtau, * R+dR,aT+daT,ycm+dycm,lam+dlam,lamK+dlamK, * iflag,ypr,ytr,fs) - f minus = calqmn(T-dT,vT-dvT,et0-det0,muB-dmuB,tau-dtau, * R-dR,aT-daT,ycm-dycm,lam-dlam,lamK-dlamK, * iflag,ypr,ytr,fs) - f if ((plus.lt.0.).and.(minus.gt.0.)) then fact = plus plus = minus minus = fact elseif ((plus.ge.0.).and.(minus.le.0.)) then else write(*,*) 'error estimates below are unreliable' write(100,*) 'error estimates below are unreliable' endif endif if (i.eq.1) then write(*,4001) f,plus,minus write(100,4001) f,plus,minus elseif (i.eq.2) then write(*,4002) f,plus,minus write(100,4002) f,plus,minus elseif (i.eq.3) then write(*,4003) f,plus,minus write(100,4003) f,plus,minus elseif (i.eq.4) then write(*,4004) f,plus,minus write(100,4004) f,plus,minus elseif (i.eq.5) then write(*,4005) f,plus,minus write(100,4005) f,plus,minus elseif (i.eq.6) then write(*,4006) f,plus,minus write(100,4006) f,plus,minus elseif (i.eq.7) then write(*,4007) f,plus,minus write(100,4007) f,plus,minus elseif (i.eq.8) then write(*,4008) f,plus,minus write(100,4008) f,plus,minus elseif (i.eq.9) then write(*,4009) f,plus,minus write(100,4009) f,plus,minus elseif (i.eq.10) then write(*,4010) f,plus,minus write(100,4010) f,plus,minus elseif (i.eq.11) then write(*,4011) f,plus,minus write(100,4011) f,plus,minus elseif (i.eq.12) then write(*,4012) f,plus,minus write(100,4012) f,plus,minus elseif (i.eq.13) then write(*,4013) f,plus,minus write(100,4013) f,plus,minus endif enddo 3999 format (' vs/c = ', f12.7, ' +',2f13.7) 4000 format (' vL/c = ', f12.7, ' +',2f13.7) 4001 format (' muB = ', f12.7, ' +',2f13.7, ' MeV') 4002 format (' t1 = ', f12.7, ' +',2f13.7, ' fm/c') 4003 format (' t3 = ', f12.7, ' +',2f13.7, ' fm/c') 4004 format ('t2-t1 = ', f12.7, ' +',2f13.7, ' fm/c') 4005 format (' dtau = ', f12.7, ' +',2f13.7, ' fm/c') 4006 format (' z3 = ', f12.7, ' +',2f13.7, ' fm') 4007 format (' muS = ', f12.7, ' +',2f13.7, ' MeV') 4008 format (' muI = ', f12.7, ' +',2f13.7, ' MeV') 4009 format (' Bar# = ', f12.7, ' +',2f13.7) 4010 format ('Bproj = ', f12.7, ' +',2f13.7) 4011 format ('Btarg = ', f12.7, ' +',2f13.7) 4012 format ('brdn1 = ', f12.8, ' +',2f13.8, ' fm^-3') 4013 format ('brdn2 = ', f12.8, ' +',2f13.8, ' fm^-3') write(*,*) write(100,*) write(*,*) 'The Covariance Matrix is:' write(100,*) 'The Covariance Matrix is:' do i=1,Npar write(*,*) (cov(i,j),j=1,Npar) write(100,*) (cov(i,j),j=1,Npar) enddo write(*,*) write(100,*) write(*,*) 'The Curvature Matrix is:' write(100,*) 'The Curvature Matrix is:' do i=1,Npar do j=1,Npar cov(i,j) = alpha(i,j) enddo write(*,*) (cov(i,j),j=1,Npar) write(100,*) (cov(i,j),j=1,Npar) enddo c endif END SUBROUTINE comp(x,param,P,dP,mparam) c c Calculates either one-particle distributions or two-particle c correlation functions along with their derivatives with respect to c the relevant parameters. c param(1) = log(T/140.) T in MeV c param(2) = .5*log(vT/(1.-vT)) c param(3) = log(et0/2.) c param(4) = log(muB) muB in MeV c param(5) = log(tau/.02027) tau in MeV^{-1} c param(6) = log(R/.03041) R in MeV^{-1} c param(7) = .5*log((1.+aT)/(1.-aT)) c param(8) = yCM c param(9) = log(fs) c param(10) = log(lambda_pi) c param(11) = log(lambda_K) c param(12) = log(normpi) c param(13) = log(normK) c param(14) = rw rw in cm c param(15) = muS muS in MeV c param(16) = muI muI in MeV c param(17) = dist_id c = 1. P = pi- distribution c = 2. P = pi+ distribution c = 3. P = K+ distribution c = 4. P = K- distribution c = 5. P = proton distribution c = 6. P = antiproton distribution c = 7. P = pi- correlation function c = 8. P = pi+ correlation function c = 9. P = K+ correlation function c = 10. P = K- correlation function c IMPLICIT NONE COMMON /int/ n1,n2 c INTEGER mparam,Ndim,Ifirst,i,ichem,n1(6),n2(6),ii,iter PARAMETER (Ndim=5) REAL x(Ndim),param(mparam),P,dP(13),y,pT,T,vT,et0,muB,muS,muI,rw REAL p0(17),dmuS(9),dmuI(9),tau,R,lam,dy,qout,qside,pi2,aT,B REAL P2,dP2(10),lamK,dPd(10),ycm,fs,normpi,normK PARAMETER (pi2=6.2831853) SAVE p0,T,vT,et0,muB,tau,R,aT,lam,lamK,muS,muI,dmuS,dmuI SAVE Ifirst,rw,iter DATA Ifirst,iter /1,0/ DATA p0 /17*0./ DATA T,vT,et0,muB,tau,R,aT,lam,lamK,fs,normpi,normK * /140.,.5,2.,1.,.02027,.03041,0.,1.,1.,1.,1.,1./ c ycm = param(8) y = x(1) - ycm pT = x(2) dy = x(3) qout = x(4) qside = x(5) c c Only re-evaluate parameters if they have changed c if (Ifirst.eq.1) then Ifirst = 0 ichem = 1 rw = param(14) endif if (p0(1).ne.param(1)) then ichem = 1 p0(1) = param(1) T = 140.*exp(param(1)) endif if (p0(2).ne.param(2)) then ichem = 1 p0(2) = param(2) vT = .5*(1. + tanh(param(2))) endif if (p0(3).ne.param(3)) then ichem = 1 p0(3) = param(3) et0 = 2.*exp(param(3)) endif if (p0(4).ne.param(4)) then ichem = 1 p0(4) = param(4) muB = param(4) endif if (p0(5).ne.param(5)) then ichem = 1 p0(5) = param(5) tau = .02027*exp(param(5)) endif if (p0(6).ne.param(6)) then ichem = 1 p0(6) = param(6) R = .03041*exp(param(6)) endif if (p0(7).ne.param(7)) then ichem = 1 p0(7) = param(7) aT = tanh(param(7)) endif if (p0(8).ne.ycm) then ichem = 1 p0(8) = ycm endif if (p0(9).ne.param(9)) then ichem = 1 p0(9) = param(9) fs = exp(param(9)) endif if (p0(10).ne.param(10)) then ichem = 1 p0(10) = param(10) lam = exp(param(10)) endif if (p0(11).ne.param(11)) then ichem = 1 p0(11) = param(11) lamK = exp(param(11)) endif if (p0(12).ne.param(12)) then ichem = 1 p0(12) = param(12) normpi = exp(param(12)) endif if (p0(13).ne.param(13)) then ichem = 1 p0(13) = param(13) normK = exp(param(13)) endif c c For each iteration, find muS, muI and their derivatives c if (ichem.eq.1) then ichem = 0 iter = iter + 1 write(*,*) write(100,*) write(*,*) write(100,*) write(*,*) 'ITERATION # ',iter write(100,*) 'ITERATION # ',iter c c Write out the current parameters being used c write(*,10) T,vT,et0,muB,tau*197.33,R*197.33, & aT,ycm,fs,lam,lamK, & normpi,normK write(100,10) T,vT,et0,muB,tau*197.33,R*197.33, & aT,ycm,fs,lam,lamK, & normpi,normK 10 format( &' T ',f11.7,/, &' vT ',f11.7,/, &' et0 ',f11.7,/, &' muB/T ',f11.7,/, &' tau ',f11.7,/, &' R ',f11.7,/, &' aT ',f11.7,/, &' ycm ',f11.7,/, &' fs ',f11.7,/, &' lampi ',f11.7,/, &' lamK ',f11.7,/, &' normpi ',f11.7,/, &' normK ',f11.7) c c Solve constraint equations for muS and muI c B = 0. CALL chem(B,T,vT,et0,muB,tau,R,aT,ycm,muS,muI,dmuS,dmuI,fs) param(15) = muS param(16) = muI write(*,*) '*** Calculating # of points, dof, chi^2, ', & 'chi^2/dof for each data set ...' write(100,*) '*** Calculating # of points, dof, chi^2, ', & 'chi^2/dof for each data set ...' endif c c Evaluate distributions or correlation functions c c Change sign of chemical potentials when doing pi+, K- or antiprotons c if (mod(param(17),2.).eq.0.) then muI = -muI muS = -muS muB = -muB do i=1,9 dmuS(i) = -dmuS(i) dmuI(i) = -dmuI(i) enddo endif c c pi- or pi+ distribution c if ((param(17).eq.1.).or.(param(17).eq.2)) then CALL distpim(P,dPd,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) P = normpi*P do i=1,9 dP(i) = normpi*dPd(i) enddo dP(10) = 0. dP(11) = 0. dP(12) = P dP(13) = 0. c c K+ or K- distribution c elseif ((param(17).eq.3.).or.(param(17).eq.4.)) then CALL distkp(P,dPd,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) P = normK*P do i=1,9 dP(i) = normK*dPd(i) enddo dP(10) = 0. dP(11) = 0. dP(12) = 0. dP(13) = P c c proton or antiproton distribution c elseif ((param(17).eq.5.).or.(param(17).eq.6.)) then CALL distpro(P,dPd,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) do i=1,9 dP(i) = dPd(i) enddo dP(10) = 0. dP(11) = 0. dP(12) = 0. dP(13) = 0. c c pi- or pi+ correlation function c elseif ((param(17).eq.7.).or.(param(17).eq.8.)) then CALL corpim(P,dPd,y,pT,dy,qout,qside,T,vT,et0,muB,tau,R,aT,lam, * muS,muI,dmuS,dmuI,rw,fs,n1,n2) if (n2(6).lt.0) then CALL corpim(P2,dP2,y,pT,dy,qout,qside,T,vT,et0,muB,tau,R, * aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) P = .5*(P+P2) do ii=1,10 dPd(ii) = .5*(dPd(ii)+dP2(ii)) enddo endif do i=1,10 dP(i) = dPd(i) enddo dP(11) = 0. dP(12) = 0. dP(13) = 0. c c K+ or K- correlation function c elseif ((param(17).eq.9.).or.(param(17).eq.10.)) then CALL corkp(P,dPd,y,pT,dy,qout,qside,T,vT,et0,muB,tau,R,aT,lamK, * muS,muI,dmuS,dmuI,rw,fs,n1,n2) if (n2(6).lt.0) then CALL corkp(P2,dP2,y,pT,-dy,qout,qside,T,vT,et0,muB,tau,R,aT, * lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) P = .5*(P+P2) do ii=1,10 dPd(ii) = .5*(dPd(ii)+dP2(ii)) enddo endif do i=1,9 dP(i) = dPd(i) enddo dP(10) = 0. dP(11) = dPd(10) dP(12) = 0. dP(13) = 0. endif c c Change back signs of chemical potentials c if (mod(param(17),2.).eq.0.) then muI = -muI muS = -muS muB = -muB do i=1,9 dmuS(i) = -dmuS(i) dmuI(i) = -dmuI(i) enddo endif return END FUNCTION chi2(pam,x,y,sig,Ndata,param,mask) c c Calculates chi^2 = sum_i((yexp_i - ytheory_i)/sig_i)^2 c IMPLICIT NONE COMMON /simplex/chi2min,param_min INTEGER Npoints, Ndim, NP,i,j,Ndata(10),k,Ifirst,Npar,ndatot PARAMETER (Npoints=1000, Ndim=5, NP=17) INTEGER mask(NP) REAL x(10,Ndim,Npoints),y(10,Npoints),sig(10,Npoints),chi2i REAL chi2,P,dP(10),x0(Ndim),pam(NP),param(NP),fact(10),dof REAL param_min(NP),chi2min SAVE Ifirst,Npar,ndatot,fact,dof DATA Ifirst /1/ DATA chi2min /-1./ c c Calculate number of parameters, degree of freedom per data set c (first time only) c if (Ifirst.eq.1) then Ifirst = 0 Npar = 0 do k=1,13 if (mask(k).eq.1) Npar = Npar + 1 enddo ndatot = 0 do k=1,10 ndatot = ndatot + Ndata(k) enddo dof = float(ndatot-Npar) do k=1,10 fact(k) = float(Ndata(k)*(ndatot-Npar))/float(ndatot) enddo endif c c Unpack the variables in pam into the full parameter vector param c i = 1 do j=1,13 if (mask(j).eq.1) then param(j) = pam(i) i = i + 1 endif enddo c c Calculate chi^2 for each data set and total chi^2 c chi2 = 0. chi2i = 0. do k=1,10 param(17) = Real(k) do i=1,Ndata(k) do j=1,Ndim x0(j) = x(k,j,i) enddo CALL comp(x0,param,P,dP,NP) chi2 = chi2 + ((y(k,i)-P)/sig(k,i))**2 enddo if (Ndata(k).ne.0) then chi2i = chi2 - chi2i write(*,*) Ndata(k),fact(k),chi2i,chi2i/fact(k) write(100,*) Ndata(k),fact(k),chi2i,chi2i/fact(k) chi2i = chi2 endif enddo c c Save the minimum value of chi^2 yet seen and its parameters c if ((chi2.lt.chi2min).or.(chi2min.eq.-1.)) then chi2min = chi2 do j=1,13 param_min(j) = param(j) enddo endif write(*,*) ' chi^2 =',chi2, * ' chi^2/dof =',chi2/dof write(100,*) ' chi^2 =',chi2, * ' chi^2/dof =',chi2/dof return END SUBROUTINE chem(B0,T0,vT0,et00,muB0,tau0,R0,aT0,ycm, * muS,muI,dmuS,dmuI,fs0) c c Determines muS (=muS/T) and muI (=muI/T) such that c total strangeness=0 and isospin/baryon = original value c (For Pb+Pb I/B_0 = .5*(82-126)/208 = -0.10577) c c Also returns: c dmuS(1) = T*dmuS/dT dmuI(1) = T*dmuI/dT c dmuS(2) = vT*dmuS/dvT dmuI(2) = vT*dmuI/dvT c dmuS(3) = et0*dmuS/det0 = 0 = dmuI(3) = et0*dmuI/det0 c dmuS(4) = dmuS/dmuB dmuI(4) = dmuI/dmuB c dmuS(5) = tau*dmuS/dtau dmuI(5) = tau*dmuI/dtau c dmuS(6) = R*dmuS/dR dmuI(6) = R*dmuI/dR c dmuS(7) = dmuS/daT dmuI(7) = dmuI/daT c dmuS(8) = dmuS/dycm dmuI(8) = dmuI/dycm c dmuS(9) = fs*dmuS/dfs dmuI(9) = fs*dmuI/dfs c IMPLICIT NONE COMMON /asym/ Zp,Ap,Ztar,Atar,yp,yt COMMON /params/ muB,T,vT,et0,tau,R,aT,IB0,B,fs c INTEGER n,ifail PARAMETER (n=2) REAL muB,T,vT,et0,IB0,mu(n),muS,muI,T0,vT0,et00,muB0 REAL dmuS(9),dmuI(9),B,tau0,tau,R0,R,aT0,aT,fs,fs0,B0,dIB0 REAL ycm,Zp,Ap,Ztar,Atar,yp,yt,shpcm,shcmt real tolx, tolf SAVE ifail,tolx,tolf DATA ifail /0/ DATA tolx,tolf /0.000001,0.000001/ c c When getting errors on final physical quantities, do not need c as much accuracy. c if (B0.eq.-1.) then tolx = .0001 tolf = .0001 endif c B = B0 T = T0 vT = vT0 et0 = et00 muB = muB0 tau = tau0 R = R0 aT = aT0 fs = fs0 c c Calculate initial isospin per baryon and its derivative c with respect to ycm. c shcmt = sinh(ycm-yt) shpcm = sinh(yp-ycm) IB0 = -.5 + (Zp*shcmt/Ap + Ztar*shpcm/Atar)/(shpcm+shcmt) dIB0 = (Zp/Ap-Ztar/Atar)*sinh(yp-yt)/(shpcm+shcmt)**2 c c Solve for muS and muI c mu(1) = 0. mu(2) = 0. CALL mnewt(20,mu,n,tolx,tolf,ifail) muS = mu(1) muI = mu(2) if (ifail.eq.1) then if (B0.eq.-1.) then ifail=0 return else write(*,*) 'muS = ',muS,' muI = ',muI write(100,*) 'muS = ',muS,' muI = ',muI stop endif endif c c Write out projectile baryons, target baryons, total baryons, c and isospin/baryon c if (B0.ne.-1.) then write(*,500) T*muS, T*muI write(100,500) T*muS, T*muI write(*,1000) B*shcmt/(shcmt+shpcm),B*shpcm/(shcmt+shpcm) write(100,1000) B*shcmt/(shcmt+shpcm),B*shpcm/(shcmt+shpcm) write(*,2000) B,IB0 write(100,2000) B,IB0 c c Calculate derivatives dmuS and dmuI c CALL dquant(dmuS,dmuI,T,vT,et0,muB,tau,R,aT,muS,muI,fs,dIB0) endif c 500 format(' muS: ',f13.7,' muI:',f13.7) 1000 format(' Projectile baryons: ',f13.7,' Target baryons:',f13.7) 2000 format(' Total baryon number:',f13.7,' Isospin/Baryon:',f13.7) B0 = B return END SUBROUTINE usrfun(mu,n,NP,fvec,fjac) c c Called by the 2-D solver mnewt c Calculates total strangeness (S) = fvec(1) c and isospin/baryon (IB) - I/B initial (IB0) = fvec(2) c Also calculates derivatives of these with respect to muS and muI c IMPLICIT NONE COMMON /params/ muB,T,vT,et0,tau,R,aT,IB0,B,fs c INTEGER n,NP REAL mu(n),fvec(NP),fjac(NP,NP),muB,T,vT,et0,tau,R,aT REAL S,dSS,dSI,IB,dIBS,dIBI,B,IB0,fs c call quant(S,dSS,dSI,IB,dIBS,dIBI,B,T,vT,et0,muB,tau,R,aT, * mu(1),mu(2),fs) fvec(1) = S fvec(2) = IB - IB0 fjac(1,1) = dSS fjac(1,2) = dSI fjac(2,1) = dIBS fjac(2,2) = dIBI return END SUBROUTINE quant(S,dSS,dSI,IB,dIBS,dIBI,B,T,vT,et0,muB,tau,R,aT, * muS,muI,fs) c c For the given parameters, this function determines the total c strangeness (S), isospin/baryon (IB) and their derivatives with c respect to muS/T (dSS,dIBS) and muI/T (dSI,dIBI). Also baryon # (B) c IMPLICIT NONE INTEGER kb REAL T,vT,et0,muB,muS,muI,mass,S,dSS,dSI,IB,dIBS,dIBI,B,I,dIS,dII REAL dBS,dBI,shS,chS,shB,chB,sh5I,ch5I,shI,chI,shBS,chBS,tau,R,aT REAL sh15I,ch15I,shB2S,chB2S,number,x,B0,S0,I0,tol,fs SAVE tol DATA tol /0.000001/ c c When getting errors on final physical quantities, do not need c as much accuracy. c if (B.eq.-1.) tol = 0.0001 c S = 0. dSS = 0. dSI = 0. I = 0. dIS = 0. dII = 0. B = 0. dBS = 0. dBI = 0. do kb=1,15 shS = sinh(kb*muS) chS = cosh(kb*muS) shB = sinh(kb*muB) chB = cosh(kb*muB) sh5I = sinh(.5*kb*muI) ch5I = cosh(.5*kb*muI) shI = 2.*sh5I*ch5I chI = 2.*ch5I*ch5I - 1. sh15I = shI*ch5I + sh5I*chI ch15I = chI*ch5I + shI*sh5I shBS = sinh(kb*(muB-muS)) chBS = cosh(kb*(muB-muS)) shB2S = sinh(kb*(muB-2.*muS)) chB2S = cosh(kb*(muB-2.*muS)) c c pions c mass = 139.6 x = number(kb,mass,T,vT,et0,tau,R,aT) I = I + 2.*shI*x dII = dII + 2.*kb*chI*x c c Kaons and Kstars c mass = 494. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) S = S + 4.*ch5I*shS*x dSS = dSS + 4.*kb*ch5I*chS*x dSI = dSI + 2.*kb*sh5I*shS*x I = I + 2.*sh5I*chS*x dIS = dIS + 2.*kb*sh5I*shS*x dII = dII + kb*ch5I*chS*x c mass = 892. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) S = S + 12.*ch5I*shS*x dSS = dSS + 12.*kb*ch5I*chS*x dSI = dSI + 6.*kb*sh5I*shS*x I = I + 6.*sh5I*chS*x dIS = dIS + 6.*kb*sh5I*shS*x dII = dII + 3.*kb*ch5I*chS*x c c rhos c mass = 770. x = number(kb,mass,T,vT,et0,tau,R,aT) I = I + 6.*shI*x dII = dII + 6.*kb*chI*x c c nucleons and nucleon stars c mass = 939. x = number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x I = I + 4.*sh5I*chB*x dII = dII + 2.*kb*ch5I*chB*x B = B + 8.*ch5I*shB*x dBI = dBI + 4.*kb*sh5I*shB*x c c mass = 1440. c x = number(kb,mass,T,vT,et0,tau,R,aT) c if (mod(kb,2).eq.0) x = -x c I = I + 4.*sh5I*chB*x c dII = dII + 2.*kb*ch5I*chB*x c B = B + 8.*ch5I*shB*x c dBI = dBI + 4.*kb*sh5I*shB*x c c mass = 1525. c x = number(kb,mass,T,vT,et0,tau,R,aT) c if (mod(kb,2).eq.0) x = -x c I = I + 12.*sh5I*chB*x c dII = dII + 6.*kb*ch5I*chB*x c B = B + 24.*ch5I*shB*x c dBI = dBI + 12.*kb*sh5I*shB*x c c Deltas c mass = 1232. x = number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x I = I + 8.*chB*(sh5I+3.*sh15I)*x dII = dII + 4.*kb*chB*(ch5I+9.*ch15I)*x B = B + 16.*shB*(ch5I+ch15I)*x dBI = dBI + 8.*kb*shB*(sh5I+3.*sh15I)*x c c Lambdas and lambda(1405) c mass = 1116. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x S = S - 4.*shBS*x dSS = dSS + 4.*kb*chBS*x B = B + 4.*shBS*x dBS = dBS - 4.*kb*chBS*x c mass = 1405. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x S = S - 4.*shBS*x dSS = dSS + 4.*kb*chBS*x B = B + 4.*shBS*x dBS = dBS - 4.*kb*chBS*x c c Sigmas and sigma(1385) c mass = 1194. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x S = S - 4.*shBS*(1.+2.*chI)*x dSS = dSS + 4.*kb*chBS*(1.+2.*chI)*x dSI = dSI - 8.*kb*shBS*shI*x I = I + 8.*chBS*shI*x dIS = dIS - 8.*kb*shBS*shI*x dII = dII + 8.*kb*chBS*chI*x B = B + 4.*shBS*(1.+2.*chI)*x dBS = dBS - 4.*kb*chBS*(1.+2.*chI)*x dBI = dBI + 8.*kb*shBS*shI*x c mass = 1385. x = fs*number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x S = S - 8.*shBS*(1.+2.*chI)*x dSS = dSS + 8.*kb*chBS*(1.+2.*chI)*x dSI = dSI - 16.*kb*shBS*shI*x I = I + 16.*chBS*shI*x dIS = dIS - 16.*kb*shBS*shI*x dII = dII + 16.*kb*chBS*chI*x B = B + 8.*shBS*(1.+2.*chI)*x dBS = dBS - 8.*kb*chBS*(1.+2.*chI)*x dBI = dBI + 16.*kb*shBS*shI*x c c Hyperons c mass = 1318. x = fs*fs*number(kb,mass,T,vT,et0,tau,R,aT) if (mod(kb,2).eq.0) x = -x S = S - 16.*shB2S*ch5I*x dSS = dSS + 32.*kb*chB2S*ch5I*x dSI = dSI - 8.*kb*shB2S*sh5I*x I = I + 4.*chB2S*sh5I*x dIS = dIS - 8.*kb*shB2S*sh5I*x dII = dII + 2.*kb*chB2S*ch5I*x B = B + 8.*shB2S*ch5I*x dBS = dBS - 16.*kb*chB2S*ch5I*x dBI = dBI + 4.*kb*shB2S*sh5I*x c if ( (kb.gt.1).and.(abs(B-B0).le.tol*abs(B)) * .and.(abs(I-I0).le.tol*abs(I)) * .and.(abs(S-S0).le.tol*abs(S)) ) goto 999 B0 = B I0 = I S0 = S enddo c 999 IB = I/B dIBS = dIS/B - I*dBS/B**2 dIBI = dII/B - I*dBI/B**2 return END SUBROUTINE dquant(dmuS,dmuI,T,vT,et0,muB,tau,R,aT,muS,muI,fs,dIB0) c c For the given parameters, this function determines the derivatives c of muS (=muS/T) and muI (=muI/T) with respect to those parameters. c Explicitly: c dmuS(1) = T*dmuS/dT dmuI(1) = T*dmuI/dT c dmuS(2) = vT*dmuS/dvT dmuI(2) = vT*dmuI/dvT c dmuS(4) = dmuS/dmuB dmuI(4) = dmuI/dmuB , etc. c IMPLICIT NONE INTEGER kb,i REAL T,vT,et0,muB,muS,muI,mass,S,dSS,dSI,IB,dIBS,dIBI,B,dIB0 REAL dmuS(9),dmuI(9),dS(9),dI(9),dB(9),dB0(9),dI0,dS0 REAL shS,chS,shB,chB,sh5I,ch5I,shI,chI,shBS,chBS,dN(7) REAL sh15I,ch15I,shB2S,chB2S,N,denom,dIB(9),tol,tau,R,aT,fs PARAMETER (tol=.000001) c CALL quant(S,dSS,dSI,IB,dIBS,dIBI,B,T,vT,et0,muB,tau,R,aT, * muS,muI,fs) c do i=1,9 dS(i) = 0. dI(i) = 0. dB(i) = 0. enddo do kb=1,15 shS = sinh(kb*muS) chS = cosh(kb*muS) shB = sinh(kb*muB) chB = cosh(kb*muB) sh5I = sinh(.5*kb*muI) ch5I = cosh(.5*kb*muI) shI = 2.*sh5I*ch5I chI = 2.*ch5I*ch5I - 1. sh15I = shI*ch5I + sh5I*chI ch15I = chI*ch5I + shI*sh5I shBS = sinh(kb*(muB-muS)) chBS = cosh(kb*(muB-muS)) shB2S = sinh(kb*(muB-2.*muS)) chB2S = cosh(kb*(muB-2.*muS)) c c pions c mass = 139.6 CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) do i=1,7 dI(i) = dI(i) + 2.*shI*dN(i) enddo c c Kaons and Kstars c c if (kb.lt.4) then mass = 494. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) do i=1,7 dS(i) = dS(i) + 4.*fs*ch5I*shS*dN(i) dI(i) = dI(i) + 2.*fs*sh5I*chS*dN(i) enddo dS(9) = dS(9) + 4.*fs*ch5I*shS*N dI(9) = dI(9) + 2.*fs*sh5I*chS*N c endif c c if (kb.lt.3) then mass = 892. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) do i=1,7 dS(i) = dS(i) + 12.*fs*ch5I*shS*dN(i) dI(i) = dI(i) + 6.*fs*sh5I*chS*dN(i) enddo dS(9) = dS(9) + 12.*fs*ch5I*shS*N dI(9) = dI(9) + 6.*fs*sh5I*chS*N c c rhos c mass = 770. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) do i=1,7 dI(i) = dI(i) + 6.*shI*dN(i) enddo c c nucleons and nucleon stars c mass = 939. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dI(i) = dI(i) + 4.*sh5I*chB*dN(i) dB(i) = dB(i) + 8.*ch5I*shB*dN(i) enddo dI(4) = dI(4) + 4.*kb*sh5I*shB*N dB(4) = dB(4) + 8.*kb*ch5I*chB*N c c mass = 1440. c CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) c if (mod(kb,2).eq.0) then c N = -N c do i=1,7 c dN(i) = -dN(i) c enddo c endif c do i=1,7 c dI(i) = dI(i) + 4.*sh5I*chB*dN(i) c dB(i) = dB(i) + 8.*ch5I*shB*dN(i) c enddo c dI(4) = dI(4) + 4.*kb*sh5I*shB*N c dB(4) = dB(4) + 8.*kb*ch5I*chB*N c c mass = 1525. c CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) c if (mod(kb,2).eq.0) then c N = -N c do i=1,7 c dN(i) = -dN(i) c enddo c endif c do i=1,7 c dI(i) = dI(i) + 12.*sh5I*chB*dN(i) c dB(i) = dB(i) + 24.*ch5I*shB*dN(i) c enddo c dI(4) = dI(4) + 12.*kb*sh5I*shB*N c dB(4) = dB(4) + 24.*kb*ch5I*chB*N c c Deltas c mass = 1232. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dI(i) = dI(i) + 8.*chB*(sh5I+3.*sh15I)*dN(i) dB(i) = dB(i) + 16.*shB*(ch5I+ch15I)*dN(i) enddo dI(4) = dI(4) + 8.*kb*shB*(sh5I+3.*sh15I)*N dB(4) = dB(4) + 16.*kb*chB*(ch5I+ch15I)*N c c Lambdas and lambda(1405) c mass = 1116. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dS(i) = dS(i) - 4.*fs*shBS*dN(i) dB(i) = dB(i) + 4.*fs*shBS*dN(i) enddo dS(4) = dS(4) - 4.*kb*fs*chBS*N dB(4) = dB(4) + 4.*kb*fs*chBS*N dS(9) = dS(9) - 4.*fs*shBS*N dB(9) = dB(9) + 4.*fs*shBS*N c mass = 1405. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dS(i) = dS(i) - 4.*fs*shBS*dN(i) dB(i) = dB(i) + 4.*fs*shBS*dN(i) enddo dS(4) = dS(4) - 4.*kb*fs*chBS*N dB(4) = dB(4) + 4.*kb*fs*chBS*N dS(9) = dS(9) - 4.*fs*shBS*N dB(9) = dB(9) + 4.*fs*shBS*N c c Sigmas and sigma(1385) c mass = 1194. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dS(i) = dS(i) - 4.*fs*shBS*(1.+2.*chI)*dN(i) dI(i) = dI(i) + 8.*fs*chBS*shI*dN(i) dB(i) = dB(i) + 4.*fs*shBS*(1.+2.*chI)*dN(i) enddo dS(4) = dS(4) - 4.*kb*fs*chBS*(1.+2.*chI)*N dI(4) = dI(4) + 8.*kb*fs*shBS*shI*N dB(4) = dB(4) + 4.*kb*fs*chBS*(1.+2.*chI)*N dS(9) = dS(9) - 4.*fs*shBS*(1.+2.*chI)*N dI(9) = dI(9) + 8.*fs*chBS*shI*N dB(9) = dB(9) + 4.*fs*shBS*(1.+2.*chI)*N c mass = 1385. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dS(i) = dS(i) - 8.*fs*shBS*(1.+2.*chI)*dN(i) dI(i) = dI(i) + 16.*fs*chBS*shI*dN(i) dB(i) = dB(i) + 8.*fs*shBS*(1.+2.*chI)*dN(i) enddo dS(4) = dS(4) - 8.*kb*fs*chBS*(1.+2.*chI)*N dI(4) = dI(4) + 16.*kb*fs*shBS*shI*N dB(4) = dB(4) + 8.*kb*fs*chBS*(1.+2.*chI)*N dS(9) = dS(9) - 8.*fs*shBS*(1.+2.*chI)*N dI(9) = dI(9) + 16.*fs*chBS*shI*N dB(9) = dB(9) + 8.*fs*shBS*(1.+2.*chI)*N c c Hyperons c mass = 1318. CALL dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) if (mod(kb,2).eq.0) then N = -N do i=1,7 dN(i) = -dN(i) enddo endif do i=1,7 dS(i) = dS(i) - 16.*fs*fs*shB2S*ch5I*dN(i) dI(i) = dI(i) + 4.*fs*fs*chB2S*sh5I*dN(i) dB(i) = dB(i) + 8.*fs*fs*shB2S*ch5I*dN(i) enddo dS(4) = dS(4) - 16.*kb*fs*fs*chB2S*ch5I*N dI(4) = dI(4) + 4.*kb*fs*fs*shB2S*sh5I*N dB(4) = dB(4) + 8.*kb*fs*fs*chB2S*ch5I*N dS(9) = dS(9) - 32.*fs*fs*shB2S*ch5I*N dI(9) = dI(9) + 8.*fs*fs*chB2S*sh5I*N dB(9) = dB(9) + 16.*fs*fs*shB2S*ch5I*N c endif c if ( (kb.gt.1).and.(abs(dI(1)-dI0).le.tol*abs(dI(1))) * .and.(abs(dS(1)-dS0).le.tol*abs(dS(1))) * .and.(abs(dB(1)-dB0(1)).le.tol*abs(dB(1))) * .and.(abs(dB(2)-dB0(2)).le.tol*abs(dB(2))) * .and.(abs(dB(4)-dB0(4)).le.tol*abs(dB(4))) * .and.(abs(dB(5)-dB0(5)).le.tol*abs(dB(5))) * .and.(abs(dB(6)-dB0(6)).le.tol*abs(dB(6))) * .and.(abs(dB(7)-dB0(7)).le.tol*abs(dB(7))) ) goto 999 if (kb.eq.15) then write(*,*) 'WARNING: full 15 terms in kb sum used in dquant' write(100,*) 'WARNING: full 15 terms in kb sum used in dquant' write(*,*) 'dI(1),dI0',dI(1),dI0,(dI(1)-dI0)/dI(1) write(100,*) 'dI(1),dI0',dI(1),dI0 write(*,*) 'dS(1),dS0',dS(1),dS0,(dS(1)-dS0)/dS(1) write(100,*) 'dS(1),dS0',dS(1),dS0 do i=1,7 write(*,998) i,i,dB(i),dB0(i),(dB(i)-dB0(i))/dB(i) write(100,998) i,i,dB(i),dB0(i),(dB(i)-dB0(i))/dB(i) enddo 998 format(' dB(',i1,'),dB0(',i1,')',3g16.7) endif dI0 = dI(1) dS0 = dS(1) do i=1,7 dB0(i) = dB(i) enddo enddo c 999 denom = dSS*dIBI - dSI*dIBS c c T*dmu/dT dIB(1) = (dI(1)+IB*B - IB*(dB(1)+B))/B dmuS(1) = (dSI*dIB(1) - dIBI*dS(1))/denom dmuI(1) = (dIBS*dS(1) - dSS*dIB(1))/denom c c vT*dmu/dvT dIB(2) = (dI(2) - IB*dB(2))/B dmuS(2) = (dSI*dIB(2) - dIBI*dS(2))/denom dmuI(2) = (dIBS*dS(2) - dSS*dIB(2))/denom c c dmu/det0 = 0 c c dmu/dmuB dIB(4) = (dI(4) - IB*dB(4))/B dmuS(4) = (dSI*dIB(4) - dIBI*dS(4))/denom dmuI(4) = (dIBS*dS(4) - dSS*dIB(4))/denom c c tau*dmu/dtau dIB(5) = (dI(5)+IB*B - IB*(dB(5)+B))/B dmuS(5) = (dSI*dIB(5) - dIBI*dS(5))/denom dmuI(5) = (dIBS*dS(5) - dSS*dIB(5))/denom c c R*dmu/dR dIB(6) = (dI(6)+2.*IB*B - IB*(dB(6)+2.*B))/B dmuS(6) = (dSI*dIB(6) - dIBI*dS(6))/denom dmuI(6) = (dIBS*dS(6) - dSS*dIB(6))/denom c c dmu/daT dIB(7) = (dI(7) - IB*dB(7))/B dmuS(7) = (dSI*dIB(7) - dIBI*dS(7))/denom dmuI(7) = (dIBS*dS(7) - dSS*dIB(7))/denom c c dmu/dycm dmuS(8) = -dSI*dIB0/denom dmuI(8) = dSS*dIB0/denom c c fs*dmu/dfs dIB(9) = (dI(9) - IB*dB(9))/B dmuS(9) = (dSI*dIB(9) - dIBI*dS(9))/denom dmuI(9) = (dIBS*dS(9) - dSS*dIB(9))/denom c return END SUBROUTINE dnumber(kb,mass,T,vT,et0,tau,R,aT,N,dN) c c sum_kb exp(kb*mu/T)*N c = (total # of particles of given species) c Also: c dN(1) = T*dN/dT dN(2) = vT*dN/dvT, etc. c c The Bose-Einstein (or Fermi-Dirac) distribution is turned into a c sum over Boltzmann modes kb. The sum is taken outside the function. c The hypersurface delta function is used to do an analytic tau c integration making tau -> tau_0. c An analytic phi integration produces I_0(kb*gamT*pT*vT*(rho/R)/T). c The derivative of I_0(x) is I_1(x) c An analytic y-eta integration produces K_1(kb*gamT*mT/T). c The derivative of K_1(x) is -(K_0(x) + K_1(x)/x). c An analytic eta integration gives the factor 2*etamax where c etamax = et0*sqrt(1-(rho/R)^2) c The numerical rho integration is done using nk-point Gauss-Jacobi c quadrature with the integration variable u = 2(rho/R)^2 - 1 c and the weight function W(u) = sqrt(1-u). c The numerical mT integration is done using nl-point Gauss-Laguerre c quadrature with the integration variable w = dxj*(mT-mass) c The coefficient dxj is found by determining the slope of the log c of the distribution between mT-mass = 1500. and 2000. c IMPLICIT NONE INTEGER ii,j,k,nl,nk,kb,Ifirst PARAMETER (nl=20,nk=20) REAL q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9,b REAL s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,d REAL xl(nl),wl(nl),xk(nk),wk(nk),N,dN(7),expr,gamT,tau,R,aT REAL mass,T,vT,et0,a,c,sb(5),sc(5),pT,sd1,sd2,sd3,sd4,xj,dxj REAL mT,vT2,rhop2,rhop,pi2,fact,c1,c2,sqr REAL besi0,besk1,bessi0,bessi1,bessk0,bessk1,besi1,besk0,root2 PARAMETER (pi2=6.2831853, root2=1.41421356) SAVE q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9 SAVE s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,xk,wk,xl,wl,Ifirst DATA Ifirst /1/ c I_0 DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228e0,0.1328592e-1, * 0.225319e-2,-0.157565e-2,0.916281e-2,-0.2057706e-1, * 0.2635537e-1,-0.1647633e-1,0.392377e-2/ c I_1 DATA r1,r2,r3,r4,r5,r6,r7,r8,r9/0.39894228,-0.3988024e-1, * -0.362018e-2,0.163801e-2,-0.1031555e-1,0.2282967e-1, * -0.2895312e-1,0.1787654e-1,-0.420059e-2/ c K_0 DATA t1,t2,t3,t4,t5,t6,t7/1.25331414e0,-0.7832358e-1,0.2189568e-1, * -0.1062446e-1,0.587872e-2,-0.251540e-2,0.53208e-3/ c K_1 DATA s1,s2,s3,s4,s5,s6,s7/1.25331414,0.23498619,-0.3655620e-1, * 0.1504268e-1,-0.780353e-2,0.325614e-2,-0.68245e-3/ c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 call gaulag(xl,wl,nl,0.) call gaujac(xk,wk,nk,.5,0.) endif c vT2 = vT*vT c c Integrate over xl(j) = dxj*(mT-m) c do ii=1,5 sb(ii) = 0. enddo do j=-1,nl c c j=-1 and j=0 determine the log of distribution at mT-m=1500 and 2000. c From these values, the slope dxj is determined. c if (j.le.0) then xj = 0. if (j.eq.-1) then mT = 1500. + mass else mT = 2000. + mass endif else xj = xl(j) mT = xj/dxj + mass endif pT = sqrt(mT**2-mass**2) c c Integrate over xk(k) = 2.*(rho/R)**2 - 1. c do ii=1,5 sc(ii) = 0. enddo do k=1,nk rhop2 = .5*(xk(k) + 1.) rhop = sqrt(rhop2) gamT = 1./sqrt(1.-vT2*rhop2) a = kb*gamT*pT*vT*rhop/T c = kb*gamT*mT/T sqr = sqrt(1. + aT*rhop2) c c Calculate I_0(a)*K_1(c), I_1(a)*K_1(c), I_0(a)*K_0(c), I_1(a)*K_0(c) c if (a.lt.3.75) then besi0 = bessi0(a) besk1 = bessk1(c) besi1 = bessi1(a) besk0 = bessk0(c) expr = exp(xj) else b = 3.75/a d = 2./c besi0 = * q1+b*(q2+b*(q3+b*(q4+b*(q5+b*(q6+b*(q7+b*(q8+b*q9))))))) besi1 = * r1+b*(r2+b*(r3+b*(r4+b*(r5+b*(r6+b*(r7+b*(r8+b*r9))))))) besk1 = s1+d*(s2+d*(s3+d*(s4+d*(s5+d*(s6+d*s7))))) besk0 = t1+d*(t2+d*(t3+d*(t4+d*(t5+d*(t6+d*t7))))) expr = exp(xj+a-c)/sqrt(a*c) endif sd1 = expr*besi0*besk1 sd2 = expr*besi1*besk1 sd3 = expr*besi0*besk0 sd4 = expr*besi1*besk0 c c Evaluate integrands for N, dNdT and dNdvT c c1 = wk(k)*sqr*mT c2 = wk(k)*tau*rhop*pT/R c N sc(1) = sc(1) + c1*sd1 - aT*c2*sd4 if (j.gt.0) then c T*dN/dT - N sc(2) = sc(2) + c1*(c*sd3 - a*sd2) - aT*c2*(c*sd2 -a*sd3) c vT*dN/dvT sc(3) = sc(3) +gamT**2*( c1*(a*sd2-vT2*rhop2*(sd1+c*sd3)) + - aT*c2*(a*sd3-sd4-c*vT2*rhop2*sd2) ) c tau*dN/dtau - N = - (R*dN/dR - 2*N) sc(4) = sc(4) - aT*c2*sd4 c dN/daT sc(5) = .5*rhop2*wk(k)*mT*sd1/sqr - c2*sd4 endif enddo if (j.eq.-1) then dxj = log(mT*sc(1)) elseif (j.eq.0) then dxj = (dxj-log(mT*sc(1)))/500. else do ii=1,5 sb(ii) = sb(ii) + wl(j)*mT*sc(ii) enddo endif enddo c fact = tau*R*R*et0/(dxj*pi2*root2) N = fact*sb(1) dN(1) = fact*sb(2) dN(2) = fact*sb(3) dN(3) = 0. dN(4) = 0. dN(5) = fact*sb(4) dN(6) = -dN(5) dN(7) = fact*sb(5) if (N.gt.1.e-307) then continue else N = 0. do ii=1,7 dN(ii) = 0. enddo endif return END FUNCTION number(kb,mass,T,vT,et0,tau,R,aT) c c sum_kb exp(kb*mu/T)*N c = (total # of particles of given species) c c The Bose-Einstein (or Fermi-Dirac) distribution is turned into a c sum over Boltzmann modes kb. The sum is taken outside the function. c The hypersurface delta function is used to do an analytic tau c integration making tau -> tau_0. c An analytic phi integration produces I_0(kb*gamT*pT*vT*(rho/R)/T). c The derivative of I_0(x) is I_1(x) c An analytic y-eta integration produces K_1(kb*gamT*mT/T). c The derivative of K_1(x) is -(K_0(x) + K_1(x)/x). c An analytic eta integration gives the factor 2*etamax where c etamax = et0*sqrt(1-(rho/R)^2) c The numerical rho integration is done using nk-point Gauss-Jacobi c quadrature with the integration variable u = 2(rho/R)^2 - 1 c and the weight function W(u) = sqrt(1-u). c The numerical mT integration is done using nl-point Gauss-Laguerre c quadrature with the integration variable w = dxj*(mT-mass) c The coefficient dxj is found by determining the slope of the log c of the distribution between mT-mass = 1500. and 2000. c IMPLICIT NONE INTEGER j,k,nl,nk,kb,Ifirst PARAMETER (nl=20,nk=20) REAL q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9,b REAL s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,d REAL xl(nl),wl(nl),xk(nk),wk(nk),expr,gamT,tau,R,aT REAL number,mass,T,vT,et0,a,c,sb,sc,pT,xj,dxj REAL mT,vT2,rhop2,rhop,pi2 REAL besi0,besk1,bessi0,bessi1,bessk0,bessk1,root2,besi1,besk0 PARAMETER (pi2=6.2831853, root2=1.41421356) SAVE q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9 SAVE s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,xk,wk,xl,wl,Ifirst DATA Ifirst /1/ c I_0 DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228e0,0.1328592e-1, * 0.225319e-2,-0.157565e-2,0.916281e-2,-0.2057706e-1, * 0.2635537e-1,-0.1647633e-1,0.392377e-2/ c I_1 DATA r1,r2,r3,r4,r5,r6,r7,r8,r9/0.39894228,-0.3988024e-1, * -0.362018e-2,0.163801e-2,-0.1031555e-1,0.2282967e-1, * -0.2895312e-1,0.1787654e-1,-0.420059e-2/ c K_0 DATA t1,t2,t3,t4,t5,t6,t7/1.25331414e0,-0.7832358e-1,0.2189568e-1, * -0.1062446e-1,0.587872e-2,-0.251540e-2,0.53208e-3/ c K_1 DATA s1,s2,s3,s4,s5,s6,s7/1.25331414,0.23498619,-0.3655620e-1, * 0.1504268e-1,-0.780353e-2,0.325614e-2,-0.68245e-3/ c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 call gaulag(xl,wl,nl,0.) call gaujac(xk,wk,nk,.5,0.) endif c vT2 = vT*vT c c Integrate over xl(j) = dxj*(mT-m) c sb = 0. do j=-1,nl c c j=-1 and j=0 determine the log of distribution at mT-m=1500 and 2000. c From these values, the slope dxj is determined. c if (j.le.0) then xj = 0. if (j.eq.-1) then mT = 1500. + mass else mT = 2000. + mass endif else xj = xl(j) mT = xj/dxj + mass endif pT = sqrt(mT**2-mass**2) c c Integrate over xk(k) = 2.*(rho/R)**2 - 1. c sc = 0. do k=1,nk rhop2 = .5*(xk(k) + 1.) rhop = sqrt(rhop2) gamT = 1./sqrt(1.-vT2*rhop2) a = kb*gamT*pT*vT*rhop/T c = kb*gamT*mT/T c c Calculate I_0(a)*K_1(c), I_1(a)*K_1(c), I_0(a)*K_0(c), I_1(a)*K_0(c) c if (a.lt.3.75) then besi0 = bessi0(a) besk1 = bessk1(c) besi1 = bessi1(a) besk0 = bessk0(c) expr = exp(xj) else b = 3.75/a d = 2./c besi0 = * q1+b*(q2+b*(q3+b*(q4+b*(q5+b*(q6+b*(q7+b*(q8+b*q9))))))) besi1 = * r1+b*(r2+b*(r3+b*(r4+b*(r5+b*(r6+b*(r7+b*(r8+b*r9))))))) besk1 = s1+d*(s2+d*(s3+d*(s4+d*(s5+d*(s6+d*s7))))) besk0 = t1+d*(t2+d*(t3+d*(t4+d*(t5+d*(t6+d*t7))))) expr = exp(xj+a-c)/sqrt(a*c) endif c c Evaluate integrand for N c sc = sc + wk(k)*expr*( mT*sqrt(1.+aT*rhop2)*besi0*besk1 + - aT*pT*rhop*tau*besi1*besk0/R ) enddo if (j.eq.-1) then dxj = log(mT*sc) elseif (j.eq.0) then dxj = (dxj-log(mT*sc))/500. else sb = sb + wl(j)*mT*sc endif enddo c number = sb*tau*R*R*et0/(dxj*pi2*root2) if (number.gt.1.e-307) then continue else number = 0. endif return END SUBROUTINE hyperr(T,vT,et0,muB,tau,R,aT,ycm,lam,lamK,U0,eigen0, * Npar0,mask0,iflag0,f,plus,minus,yp0,yt0,fs) c c This subroutine uses the multidimensional solver amoeba2 to find c the maximum and minimum values obtained by certain quantities of c interest on the confidence hyperellipse. iflag0 specifies which c quantities to determine in the function calqmn. c IMPLICIT NONE COMMON /erreal/ var,U,eigen,theta,ypr,ytr COMMON /erint/ Npar,mask,iflag,isign c INTEGER NP,NP2,Npar,iflag,iflag0,i,j,iter,isign,Npar0 PARAMETER (NP=17,NP2=NP**2) INTEGER mask(NP),mask0(NP),Ifirst REAL eigen(NP),vcalq(NP),eigen0(NP),ypr,ytr,yp0,yt0,fs REAL T,vT,et0,muB,tau,R,aT,ycm,lam,lamK,U(NP,NP),U0(NP,NP) REAL var(NP),ftol,hpi,phi(NP-1),vertex(NP,NP-1) REAL calq,vcalqmn,theta(NP),plus,minus,calqmn,f EXTERNAL calq PARAMETER(ftol=0.0001,hpi=1.570796327) SAVE Ifirst DATA Ifirst /1/ c iflag = iflag0 c c Read in the eigenvalues and eigenvectors of the covariance matrix c at the minimum of chi^2. Also read in the parameters corresponding c to that minimum. c if (Ifirst.eq.1) then Ifirst = 0 Npar = Npar0 do i=1,13 mask(i) = mask0(i) enddo ypr = yp0 ytr = yt0 do i=1,Npar eigen(i) = eigen0(i) do j=1,Npar U(i,j) = U0(i,j) enddo enddo var(1) = T var(2) = vT var(3) = et0 var(4) = muB var(5) = tau var(6) = R var(7) = aT var(8) = ycm var(9) = fs var(10) = lam var(11) = lamK endif c c Evaluate "calqmn" at the minimum of chi^2 c f = calqmn(T,vT,et0,muB,tau,R,aT,ycm,lam,lamK, * iflag,ypr,ytr,fs) if (Npar.le.1) return c c Zero out the vector phi and matrix vertex c do j=1,Npar-1 phi(j) = 0. do i=1,Npar vertex(i,j) = 0. enddo enddo c c The input simplex for minimization is found by evaluating "calq" c at the positive intersection of each variable axis with the hyperellipse. c isign =1 means to minimize. c isign = 1 do i=1,Npar do j=1,i-1 vertex(i,j) = hpi phi(j) = hpi enddo vcalq(i) = calq(phi) enddo c c Minimize the function "calq" on the confidence hyperellipse c CALL amoeba2(vertex,vcalq,NP,NP-1,Npar-1,ftol,calq,iter) c c Find the minimum of the final simplex points c vcalqmn = vcalq(1) do i=2,Npar if (vcalq(i).lt.vcalqmn) then vcalqmn = vcalq(i) endif enddo minus = vcalqmn - f c c Zero out the vector phi and matrix vertex again c do j=1,Npar-1 phi(j) = 0. do i=1,Npar vertex(i,j) = 0. enddo enddo c c The input simplex for maximization is again found by evaluating "calq" c at the positive intersection of each variable axis with the hyperellipse. c isign = -1 means to maximize. c isign = -1 do i=1,Npar do j=1,i-1 vertex(i,j) = hpi phi(j) = hpi enddo vcalq(i) = calq(phi) enddo c c Maximize the function "calq" on the confidence hyperellipse c CALL amoeba2(vertex,vcalq,NP,NP-1,Npar-1,ftol,calq,iter) c c Find the maximum of the final simplex points c (actually the minimum of the negatives) c vcalqmn = vcalq(1) do i=2,Npar if (vcalq(i).lt.vcalqmn) then vcalqmn = vcalq(i) endif enddo plus = -vcalqmn - f END FUNCTION calq(phi) c c Transforms the minimization variables phi which parametrize the c surface of the hyperellipse into the standard parameters T,vT,etc. c Then calls calqmn to find the value of the desired function for c those parameters. c IMPLICIT NONE COMMON /erreal/ var,U,eigen,theta,ypr,ytr COMMON /erint/ Npar,mask,iflag,isign c INTEGER Npar,iflag,IFirst,NP,i,j,isign,isignsav PARAMETER (NP=17) INTEGER mask(NP) REAL var(NP),U(NP,NP),eigen(NP),T,vT,et0,muB,tau,R,aT,ycm,lam REAL lamK,phi(NP-1),thetil(NP),theta(NP),calq,calqmn REAL ypr,ytr,fs,calqmin SAVE Ifirst,isignsav,T,vT,et0,muB,tau,R,aT,ycm,fs,lam,lamK DATA isignsav /-2/ DATA Ifirst /1/ if (Ifirst.eq.1) then Ifirst = 0 c c These are the values of the parameters at the minimum of chi^2 c T = var(1) vT = var(2) et0 = var(3) muB = var(4) tau = var(5) R = var(6) aT = var(7) ycm = var(8) fs = var(9) lam = var(10) lamK = var(11) endif c c Convert from angular parameters to eigenvector parameters c do i=1,Npar-1 thetil(i) = cos(phi(i))/eigen(i) do j=1,i-1 thetil(i) = thetil(i)*sin(phi(j)) enddo enddo thetil(Npar) = 1./eigen(Npar) do j=1,Npar-1 thetil(Npar) = thetil(Npar)*sin(phi(j)) enddo c c Transform back to variations of normal parameters c do i=1,Npar theta(i) = 0. do j=1,Npar theta(i) = theta(i) + U(i,j)*thetil(j) enddo enddo c c Add variations theta(i) to values at the chi^2 minimum var(i) c to get trial parameter values c i = 0 if (mask(1).eq.1) then i = i + 1 T = var(1) + theta(i) if (T.le.0.) T = 0.0001 endif if (mask(2).eq.1) then i = i + 1 vT = var(2) + theta(i) if (vT.le.0.) vT = 0.0001 if (vT.ge.1.) vT = 0.9999 endif if (mask(3).eq.1) then i = i + 1 et0 = var(3) + theta(i) if (et0.le.0.) et0 = 0.0001 endif if (mask(4).eq.1) then i = i + 1 muB = var(4) + theta(i) endif if (mask(5).eq.1) then i = i + 1 tau = var(5) + theta(i) if (tau.le.0.) tau = 0.0001 endif if (mask(6).eq.1) then i = i + 1 R = var(6) + theta(i) if (R.le.0.) R = 0.0001 endif if (mask(7).eq.1) then i = i + 1 aT = var(7) + theta(i) if (aT.lt.-0.99999) aT = -0.9999 if (aT.gt.0.99999) aT = 0.9999 endif if (mask(8).eq.1) then i = i + 1 ycm = var(8) + theta(i) endif if (mask(9).eq.1) then i = i + 1 fs = var(9) + theta(i) if (fs.le.0.) fs = 0.0001 endif if (mask(10).eq.1) then i = i + 1 lam = var(10) + theta(i) if (lam.le.0.) lam = 0.0001 endif if (mask(11).eq.1) then i = i + 1 lamK = var(11) + theta(i) if (lamK.le.0.) lamK = 0.0001 endif c c Calculation of the desired quantity c calq = calqmn(T,vT,et0,muB,tau,R,aT,ycm,lam,lamK, * iflag,ypr,ytr,fs) c c If isign=-1, we maximize by minimizing the negative of calq c if (isign.eq.-1) calq = -calq c c On the first time, calqmin must be calq c if (isignsav.ne.isign) then calqmin = calq isignsav = isign endif if (calq.lt.calqmin) then calqmin = calq endif return END FUNCTION calqmn(T,vT,et0,muB,tau,R,aT,ycm,lam,lamK, * iflag,ypr,ytr,fs) c c Calculates the function to minimize on the confidence hyperellipse c The function to minimize is determined by iflag. c IMPLICIT NONE INTEGER iflag REAL T,vT,et0,muB,tau,R,aT,ycm,lam,lamK,calqmn REAL B,muS,muI,dmuS(9),dmuI(9),ypr,ytr,fs,fac,bardns PARAMETER(fac=197.33) c c Calculation of the desired quantity c c Diagnostic test cases c if (iflag.eq.-1) then calqmn = -1 elseif (iflag.eq.-2) then calqmn = vT elseif (iflag.eq.-4) then calqmn = muB c c Real cases c elseif (iflag.eq.1) then calqmn = muB*T elseif (iflag.eq.2) then calqmn = tau*sqrt(1.+aT) elseif (iflag.eq.3) then calqmn = tau*cosh(et0) elseif (iflag.eq.4) then calqmn = tau*abs(1.-sqrt(1.+aT)) elseif (iflag.eq.5) then calqmn = tau*(1.-sqrt((1.+aT)*(1.-vT*vT))) elseif (iflag.eq.6) then calqmn = tau*sinh(et0) elseif (iflag.gt.6) then c c B=-1 is a flag to tell chem not to print anything out. c B = -1. CALL chem(B,T,vT,et0,muB,tau/fac,R/fac,aT,ycm, * muS,muI,dmuS,dmuI,fs) if (iflag.eq.7) then calqmn = muS*T elseif (iflag.eq.8) then calqmn = muI*T elseif (iflag.eq.9) then calqmn = B elseif (iflag.eq.10) then calqmn = B/(1.+sinh(ypr-ycm)/sinh(ycm-ytr)) elseif (iflag.eq.11) then calqmn = B/(1.+sinh(ycm-ytr)/sinh(ypr-ycm)) elseif (iflag.eq.12) then calqmn = bardns(T,vT,muB,tau,R,aT,1.,muS,muI,fs) elseif (iflag.eq.13) then calqmn = bardns(T,vT,muB,tau,R,aT,0.,muS,muI,fs) endif endif return END FUNCTION densty(kb,mass,T,vT,tau,R,aT,rhop) c c sum_kb exp(kb*mu/T)*N c = (local density of particles of given species) c c The Bose-Einstein (or Fermi-Dirac) distribution is turned into a c sum over Boltzmann modes kb. The sum is taken outside the function. c The hypersurface delta function is used to do an analytic tau c integration making tau -> tau_0. c An analytic phi_p integration produces I_0(kb*gamT*pT*vT*(rho/R)/T). c An analytic rapidity integration produces K_1(kb*gamT*mT/T). c The numerical mT integration is done using nl-point Gauss-Laguerre c quadrature with the integration variable w = dxj*(mT-mass) c The coefficient dxj is found by determining the slope of the log c of the distribution between mT-mass = 1500. and 2000. c IMPLICIT NONE INTEGER j,nl,kb,Ifirst PARAMETER (nl=20) REAL q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9,b REAL s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,d REAL xl(nl),wl(nl),expr,gamT,tau,R,aT REAL densty,mass,T,vT,a,c,sb,sc,pT,xj,dxj REAL mT,vT2,rhop2,rhop,pi2,fac REAL besi0,besk1,bessi0,bessi1,bessk0,bessk1,besi1,besk0 PARAMETER (pi2=6.2831853, fac=197.33) SAVE q1,q2,q3,q4,q5,q6,q7,q8,q9,r1,r2,r3,r4,r5,r6,r7,r8,r9 SAVE s1,s2,s3,s4,s5,s6,s7,t1,t2,t3,t4,t5,t6,t7,xl,wl,Ifirst DATA Ifirst /1/ c I_0 DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228e0,0.1328592e-1, * 0.225319e-2,-0.157565e-2,0.916281e-2,-0.2057706e-1, * 0.2635537e-1,-0.1647633e-1,0.392377e-2/ c I_1 DATA r1,r2,r3,r4,r5,r6,r7,r8,r9/0.39894228,-0.3988024e-1, * -0.362018e-2,0.163801e-2,-0.1031555e-1,0.2282967e-1, * -0.2895312e-1,0.1787654e-1,-0.420059e-2/ c K_0 DATA t1,t2,t3,t4,t5,t6,t7/1.25331414e0,-0.7832358e-1,0.2189568e-1, * -0.1062446e-1,0.587872e-2,-0.251540e-2,0.53208e-3/ c K_1 DATA s1,s2,s3,s4,s5,s6,s7/1.25331414,0.23498619,-0.3655620e-1, * 0.1504268e-1,-0.780353e-2,0.325614e-2,-0.68245e-3/ c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 call gaulag(xl,wl,nl,0.) endif c vT2 = vT*vT rhop2 = rhop*rhop gamT = 1./sqrt(1.-vT2*rhop2) c c Integrate over xl(j) = dxj*(mT-m) c sb = 0. do j=-1,nl c c j=-1 and j=0 determine the log of distribution at mT-m=1500 and 2000. c From these values, the slope dxj is determined. c if (j.le.0) then xj = 0. if (j.eq.-1) then mT = 1500. + mass else mT = 2000. + mass endif else xj = xl(j) mT = xj/dxj + mass endif pT = sqrt(mT**2-mass**2) a = kb*gamT*pT*vT*rhop/T c = kb*gamT*mT/T c c Calculate I_0(a)*K_1(c), I_1(a)*K_1(c), I_0(a)*K_0(c), I_1(a)*K_0(c) c if (a.lt.3.75) then besi0 = bessi0(a) besk1 = bessk1(c) besi1 = bessi1(a) besk0 = bessk0(c) expr = exp(xj) else b = 3.75/a d = 2./c besi0 = * q1+b*(q2+b*(q3+b*(q4+b*(q5+b*(q6+b*(q7+b*(q8+b*q9))))))) besi1 = * r1+b*(r2+b*(r3+b*(r4+b*(r5+b*(r6+b*(r7+b*(r8+b*r9))))))) besk1 = s1+d*(s2+d*(s3+d*(s4+d*(s5+d*(s6+d*s7))))) besk0 = t1+d*(t2+d*(t3+d*(t4+d*(t5+d*(t6+d*t7))))) expr = exp(xj+a-c)/sqrt(a*c) endif c c Evaluate integrand for N c sc = expr*( mT*besi0*besk1 + - aT*pT*rhop*tau*besi1*besk0/(R*sqrt(1.+aT*rhop2)) ) if (j.eq.-1) then dxj = log(mT*sc) elseif (j.eq.0) then dxj = (dxj-log(mT*sc))/500. else sb = sb + wl(j)*mT*sc endif enddo c densty = 2.*sb/(gamT*dxj*pi2*pi2*fac**3) return END FUNCTION bardns(T,vT,muB,tau,R,aT,rhop,muS,muI,fs) c c For the given parameters, this function determines the total c local baryon density at some point in rho/R (baryon density c does not vary with eta for this boost-invariant model. c IMPLICIT NONE INTEGER kb,iprint REAL T,vT,muB,muS,muI,mass,bardns,rhop,bardns0 REAL shB,ch5I,chI,shBS,tau,R,aT REAL ch15I,shB2S,densty,x,fs,tol PARAMETER (tol=.0001) SAVE iprint DATA iprint /0/ c bardns = 0. do kb=1,6 shB = sinh(kb*muB) ch5I = cosh(.5*kb*muI) chI = 2.*ch5I*ch5I - 1. ch15I = cosh(1.5*kb*muI) shBS = sinh(kb*(muB-muS)) shB2S = sinh(kb*(muB-2.*muS)) c c nucleons c mass = 939. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 8.*ch5I*shB*x c c c Deltas c mass = 1232. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 16.*shB*(ch5I+ch15I)*x c c Lambdas and lambda(1405) c mass = 1116. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 4.*fs*shBS*x c mass = 1405. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 4.*fs*shBS*x c c Sigmas and sigma(1385) c mass = 1194. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 4.*fs*shBS*(1.+2.*chI)*x c mass = 1385. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 8.*fs*shBS*(1.+2.*chI)*x c c Hyperons c mass = 1318. x = densty(kb,mass,T,vT,tau,R,aT,rhop) if (mod(kb,2).eq.0) x = -x bardns = bardns + 8.*shB2S*ch5I*fs*fs*x if ((kb.gt.1).and.(abs(bardns-bardns0).le.tol*abs(bardns))) * goto 999 bardns0 = bardns enddo if (iprint.le.10) then write(*,*) 'Possibly not enough terms in bardns series', * abs(bardns-bardns0)/abs(bardns) write(100,*) 'Possibly not enough terms in bardns series', * abs(bardns-bardns0)/abs(bardns) if (iprint.eq.10) then write(*,*) ' Discontinuing bardns warning message' write(100,*) ' Discontinuing bardns warning message' endif iprint = iprint+1 endif c 999 continue return END SUBROUTINE mnewt(ntrial,x,n,tolx,tolf,ifail) INTEGER n,ntrial,NP,ifail REAL tolf,tolx,x(n) PARAMETER (NP=2) CU USES lubksb,ludcmp,usrfun INTEGER i,k,indx(NP) REAL d,errf,errx,fjac(NP,NP),fvec(NP),p(NP) do 14 k=1,ntrial call usrfun(x,n,NP,fvec,fjac) errf=0. do 11 i=1,n errf=errf+abs(fvec(i)) 11 continue if(errf.le.tolf)return do 12 i=1,n p(i)=-fvec(i) 12 continue call ludcmp(fjac,n,NP,indx,d,ifail) if (ifail.eq.1) return call lubksb(fjac,n,NP,indx,p) errx=0. do 13 i=1,n errx=errx+abs(p(i)) x(i)=x(i)+p(i) 13 continue if(errx.le.tolx)return 14 continue return END SUBROUTINE lubksb(a,n,np,indx,b) INTEGER n,np,indx(n) REAL a(np,np),b(n) INTEGER i,ii,j,ll REAL sum ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0)then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return END SUBROUTINE ludcmp(a,n,np,indx,d,ifail) INTEGER n,np,indx(n),NMAX,ifail,iprint REAL d,a(np,np),TINY PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER i,imax,j,k REAL aamax,dum,sum,vv(NMAX) SAVE iprint DATA iprint /1/ d=1. do 12 i=1,n aamax=0. do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.) then if (iprint.le.10) then write(*,*) 'singular matrix in ludcmp' write(100,*) 'singular matrix in ludcmp' if (iprint.eq.10) then write(*,*) ' Discontinuing ludcmp error message' write(100,*) ' Discontinuing ludcmp error message' endif iprint = iprint + 1 endif ifail = 1 return endif vv(i)=1./aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0. do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.)a(j,j)=TINY if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return END FUNCTION bessk0(x) c Returns the modified Bessel function K_0(x) for positive real x c Uses bessio(x) REAL bessk0,x REAL bessi0 REAL p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566,0.42278420,0.23069756, * 0.3488590e-1,0.262698e-2,0.10750e-3,0.74e-5/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414,-0.7832358e-1,0.2189568e-1, * -0.1062446e-1,0.587872e-2,-0.251540e-2,0.53208e-3/ c Polynomial Fit if (x.le.2.0) then y = x*x/4.0 bessk0 = (-log(x/2.0)*bessi0(x)) * + (p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else y = (2.0/x) bessk0 = (exp(-x)/sqrt(x))* * (q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7)))))) endif return END FUNCTION bessk1(x) c Returns the modified Bessel function K_1(x) for positive real x c Uses bessi1(x) REAL bessk1,x REAL bessi1 REAL p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/1.0,0.15443144,-0.67278579, * -0.18156897,-0.1919402e-1,-0.110404e-2,-0.4686e-4/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414,0.23498619,-0.3655620e-1, * 0.1504268e-1,-0.780353e-2,0.325614e-2,-0.68245e-3/ c Polynomial Fit if (x.le.2.0) then y = x*x/4.0 bessk1 = (log(x/2.0)*bessi1(x)) * + (1.0/x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else y = (2.0/x) bessk1 = (exp(-x)/sqrt(x))* * (q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7)))))) endif return END FUNCTION bessi0(x) REAL bessi0,x REAL ax REAL p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0,3.5156229,3.0899424, *1.2067492,0.2659732,0.360768e-1,0.45813e-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228,0.1328592e-1, *0.225319e-2,-0.157565e-2,0.916281e-2,-0.2057706e-1,0.2635537e-1, *-0.1647633e-1,0.392377e-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) endif return END FUNCTION bessi1(x) c Returns the modified Bessel function I_1(x) for any real x REAL bessi1,x REAL ax REAL p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5e0,0.87890594e0,0.51498869e0, * 0.15084934e0,0.2658733e-1,0.301532e-2,0.32411e-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228e0,-0.3988024e-1, * -0.362018e-2,0.163801e-2,-0.1031555e-1,0.2282967e-1, * -0.2895312e-1,0.1787654e-1,-0.420059e-2/ c Polynomial Fit if (abs(x).le.3.75) then y = (x/3.75)**2 bessi1 = x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax = abs(x) y = 3.75/ax bessi1 = (exp(ax)/sqrt(ax))* * (q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) endif return END SUBROUTINE mrqmin(x,y,sig,ndata,a,ia,ma,covar,alpha,nca,chisq, *funcs,alamda,ycalc,chisqi) IMPLICIT NONE INTEGER ma,nca,ndata(10),ia(ma),MMAX,Ndim,Np REAL alamda,chisq,funcs,a(ma),alpha(nca,nca),covar(nca,nca) PARAMETER (MMAX=20,Np=1000,Ndim=5) REAL sig(10,Np),x(10,Ndim,Np),y(10,Np),ycalc(10,Np),chisqi(10,Np) CU USES covsrt,gaussj,mrqcof INTEGER j,k,l,m,mfit,Ifirst,ndof,kk REAL ochisq,atry(MMAX),beta(MMAX),da(MMAX) SAVE ochisq,atry,beta,da,mfit,Ifirst,ndof DATA Ifirst /1/ if (Ifirst.eq.1) then Ifirst = 0 ndof = 0 do kk=1,10 ndof = ndof + ndata(kk) enddo do kk=1,ma if (ia(kk).eq.1) ndof = ndof - 1 enddo endif if(alamda.lt.0.)then mfit=0 do 11 j=1,ma if (ia(j).ne.0) mfit=mfit+1 11 continue if (alamda.eq.-1.) alamda=0.001 call mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nca,chisq, * alamda,funcs,ycalc,chisqi) write(*,*) ' chi^2 =',chisq, * ' chi^2/dof =',chisq/float(ndof) write(100,*) ' chi^2 =',chisq, * ' chi^2/dof =',chisq/float(ndof) ochisq=chisq do 12 j=1,ma atry(j)=a(j) 12 continue endif c Scott's addition for getting curvature matrix on 1st run c if (alamda.eq.-2.) alamda = 0. j=0 do 14 l=1,ma if(ia(l).ne.0) then j=j+1 k=0 do 13 m=1,ma if(ia(m).ne.0) then k=k+1 covar(j,k)=alpha(j,k) endif 13 continue covar(j,j)=alpha(j,j)*(1.+alamda) da(j)=beta(j) endif 14 continue call gaussj(covar,mfit,nca,da,1,1) if(alamda.eq.0.)then call covsrt(covar,nca,ma,ia,mfit) return endif j=0 do 15 l=1,ma if(ia(l).ne.0) then j=j+1 atry(l)=a(l)+da(j) endif 15 continue call mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,nca,chisq, * alamda,funcs,ycalc,chisqi) write(*,*) ' chi^2 =',chisq, * ' chi^2/dof =',chisq/float(ndof) write(100,*) ' chi^2 =',chisq, * ' chi^2/dof =',chisq/float(ndof) if(chisq.lt.ochisq)then alamda=0.1*alamda ochisq=chisq j=0 do 17 l=1,ma if(ia(l).ne.0) then j=j+1 k=0 do 16 m=1,ma if(ia(m).ne.0) then k=k+1 alpha(j,k)=covar(j,k) endif 16 continue beta(j)=da(j) a(l)=atry(l) endif 17 continue else alamda=10.*alamda chisq=ochisq endif return END SUBROUTINE mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nalp,chisq, *alamda,funcs,ycalc,chisqi) IMPLICIT NONE INTEGER ma,nalp,ndata(10),ia(ma),MMAX,Ndim,jj,Np,kk PARAMETER (MMAX=20,Np=1000,Ndim=5) REAL chisq,a(ma),alpha(nalp,nalp),beta(ma),x0(Ndim) REAL sig(10,Np),x(10,Ndim,Np),y(10,Np),ycalc(10,Np),chisqi(10,Np) EXTERNAL funcs INTEGER mfit,i,j,k,l,m,Ifirst,ndatot,Npar REAL dy,sig2i,wt,ymod,dyda(MMAX),alamda,chi2i,fact(10) SAVE Ifirst,ndatot,Npar,fact DATA Ifirst /1/ if (Ifirst.eq.1) then Ifirst = 0 Npar = 0 do kk=1,ma if (ia(kk).eq.1) Npar = Npar + 1 enddo ndatot = 0 do kk=1,10 ndatot = ndatot + ndata(kk) enddo do kk=1,10 fact(kk) = float(ndata(kk)*(ndatot-Npar))/float(ndatot) enddo endif mfit=0 do 11 j=1,ma if (ia(j).ne.0) mfit=mfit+1 11 continue do 13 j=1,mfit do 12 k=1,j alpha(j,k)=0. 12 continue beta(j)=0. 13 continue chisq=0. chi2i = 0. do kk=1,10 a(17) = Real(kk) do 16 i=1,ndata(kk) do jj=1,Ndim x0(jj)=x(kk,jj,i) enddo call funcs(x0,a,ymod,dyda,ma) ycalc(kk,i) = ymod sig2i=1./(sig(kk,i)*sig(kk,i)) dy=y(kk,i)-ymod j=0 do 15 l=1,ma if(ia(l).ne.0) then j=j+1 wt=dyda(l)*sig2i k=0 do 14 m=1,l if(ia(m).ne.0) then k=k+1 alpha(j,k)=alpha(j,k)+wt*dyda(m) endif 14 continue beta(j)=beta(j)+dy*wt endif 15 continue chisq=chisq+dy*dy*sig2i chisqi(kk,i) = dy*dy*sig2i 16 continue if (ndata(kk).ne.0) then c if ((alamda.eq.-2).and.(ndata(kk).ne.0)) then chi2i = chisq - chi2i write(*,*) ndata(kk),fact(kk),chi2i,chi2i/fact(kk) write(100,*) ndata(kk),fact(kk),chi2i,chi2i/fact(kk) chi2i = chisq endif enddo do 18 j=2,mfit do 17 k=1,j-1 alpha(k,j)=alpha(j,k) 17 continue 18 continue return END SUBROUTINE gaussj(a,n,np,b,m,mp) INTEGER m,mp,n,np,NMAX REAL a(np,np),b(np,mp) REAL*4 ass(13,13) PARAMETER (NMAX=50) INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX) REAL big,dum,pivinv do i=1,n do j=1,n ass(i,j) = a(i,j) enddo enddo do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0. do 13 j=1,n if(ipiv(j).ne.1)then do 12 k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big)then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then write(*,*) 'singular matrix in gaussj' write(100,*) 'singular matrix in gaussj' do l=1,n write(*,*) (ass(l,ll),ll=1,n) write(100,*) (ass(l,ll),ll=1,n) enddo stop endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue do 15 l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum 15 continue endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.) then write(*,*) 'singular matrix in gaussj' write(100,*) 'singular matrix in gaussj' do l=1,n write(*,*) (ass(l,ll),ll=1,n) write(100,*) (ass(l,ll),ll=1,n) enddo stop endif pivinv=1./a(icol,icol) a(icol,icol)=1. do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue do 17 l=1,m b(icol,l)=b(icol,l)*pivinv 17 continue do 21 ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0. do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue do 19 l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum 19 continue endif 21 continue 22 continue do 24 l=n,1,-1 if(indxr(l).ne.indxc(l))then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return END SUBROUTINE covsrt(covar,npc,ma,ia,mfit) INTEGER ma,mfit,npc,ia(ma) REAL covar(npc,npc) INTEGER i,j,k REAL swap do 12 i=mfit+1,ma do 11 j=1,i covar(i,j)=0. covar(j,i)=0. 11 continue 12 continue k=mfit do 15 j=ma,1,-1 if(ia(j).ne.0)then do 13 i=1,ma swap=covar(i,k) covar(i,k)=covar(i,j) covar(i,j)=swap 13 continue do 14 i=1,ma swap=covar(k,i) covar(k,i)=covar(j,i) covar(j,i)=swap 14 continue k=k-1 endif 15 continue return END SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter * ,x,y0,sig,Ndata,param,mask) IMPLICIT NONE INTEGER iter,mp,ndim,np,NMAX,ITMAX REAL ftol,p(mp,np),y(mp),funk PARAMETER (NMAX=20) INTEGER Npoi, Ndi, Nparam, Ndata(10), mask(17) PARAMETER (Npoi=1000, Ndi=5, Nparam=17) REAL x(10,Ndi,Npoi),y0(10,Npoi),sig(10,Npoi),param(Nparam) EXTERNAL funk CU USES amotry,funk INTEGER i,ihi,ilo,inhi,j,m,n REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry if (iter.eq.0) then ITMAX=500 else ITMAX=iter iter=0 endif 1 do 12 n=1,ndim sum=0. do 11 m=1,ndim+1 sum=sum+p(m,n) 11 continue psum(n)=sum 12 continue 2 ilo=1 if (y(1).gt.y(2)) then ihi=1 inhi=2 else ihi=2 inhi=1 endif do 13 i=1,ndim+1 if(y(i).le.y(ilo)) ilo=i if(y(i).gt.y(ihi)) then inhi=ihi ihi=i else if(y(i).gt.y(inhi)) then if(i.ne.ihi) inhi=i endif 13 continue rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) if ((rtol.lt.ftol).or.(iter.ge.ITMAX)) then swap=y(1) y(1)=y(ilo) y(ilo)=swap do 14 n=1,ndim swap=p(1,n) p(1,n)=p(ilo,n) p(ilo,n)=swap 14 continue return endif iter=iter+2 ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0 * ,x,y0,sig,Ndata,param,mask) if (ytry.le.y(ilo)) then ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0 * ,x,y0,sig,Ndata,param,mask) else if (ytry.ge.y(inhi)) then ysave=y(ihi) ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5 * ,x,y0,sig,Ndata,param,mask) if (ytry.ge.ysave) then do 16 i=1,ndim+1 if(i.ne.ilo)then do 15 j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j)) p(i,j)=psum(j) 15 continue y(i)=funk(psum,x,y0,sig,Ndata,param,mask) endif 16 continue iter=iter+ndim goto 1 endif else iter=iter-1 endif goto 2 END FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac * ,x,y0,sig,Ndata,param,mask) IMPLICIT NONE INTEGER ihi,mp,ndim,np,NMAX REAL amotry,fac,p(mp,np),psum(np),y(mp),funk PARAMETER (NMAX=20) INTEGER Npoi, Ndi, Nparam, Ndata(10),mask(17) PARAMETER (Npoi=1000, Ndi=5, Nparam=17) REAL x(10,Ndi,Npoi),y0(10,Npoi),sig(10,Npoi),param(Nparam) EXTERNAL funk CU USES funk INTEGER j REAL fac1,fac2,ytry,ptry(NMAX) fac1=(1.-fac)/ndim fac2=fac1-fac do 11 j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 11 continue ytry=funk(ptry,x,y0,sig,Ndata,param,mask) if (ytry.lt.y(ihi)) then y(ihi)=ytry do 12 j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) 12 continue endif amotry=ytry return END SUBROUTINE amoeba2(p,y,mp,np,ndim,ftol,funk,iter) INTEGER iter,mp,ndim,np,NMAX,ITMAX REAL ftol,p(mp,np),y(mp),funk PARAMETER (NMAX=20,ITMAX=2000) EXTERNAL funk CU USES amotry2,funk INTEGER i,ihi,ilo,inhi,j,m,n REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry2 iter=0 1 do 12 n=1,ndim sum=0. do 11 m=1,ndim+1 sum=sum+p(m,n) 11 continue psum(n)=sum 12 continue 2 ilo=1 if (y(1).gt.y(2)) then ihi=1 inhi=2 else ihi=2 inhi=1 endif do 13 i=1,ndim+1 if(y(i).le.y(ilo)) ilo=i if(y(i).gt.y(ihi)) then inhi=ihi ihi=i else if(y(i).gt.y(inhi)) then if(i.ne.ihi) inhi=i endif 13 continue rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) if (rtol.lt.ftol) then swap=y(1) y(1)=y(ilo) y(ilo)=swap do 14 n=1,ndim swap=p(1,n) p(1,n)=p(ilo,n) p(ilo,n)=swap 14 continue return endif if (iter.ge.ITMAX) then write(*,*) 'WARNING: ITMAX exceeded in amoeba2' write(100,*) 'WARNING: ITMAX exceeded in amoeba2' return endif iter=iter+2 ytry=amotry2(p,y,psum,mp,np,ndim,funk,ihi,-1.0) if (ytry.le.y(ilo)) then ytry=amotry2(p,y,psum,mp,np,ndim,funk,ihi,2.0) else if (ytry.ge.y(inhi)) then ysave=y(ihi) ytry=amotry2(p,y,psum,mp,np,ndim,funk,ihi,0.5) if (ytry.ge.ysave) then do 16 i=1,ndim+1 if(i.ne.ilo)then do 15 j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j)) p(i,j)=psum(j) 15 continue y(i)=funk(psum) endif 16 continue iter=iter+ndim goto 1 endif else iter=iter-1 endif goto 2 END FUNCTION amotry2(p,y,psum,mp,np,ndim,funk,ihi,fac) INTEGER ihi,mp,ndim,np,NMAX REAL amotry2,fac,p(mp,np),psum(np),y(mp),funk PARAMETER (NMAX=20) EXTERNAL funk CU USES funk INTEGER j REAL fac1,fac2,ytry,ptry(NMAX) fac1=(1.-fac)/ndim fac2=fac1-fac do 11 j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 11 continue ytry=funk(ptry) if (ytry.lt.y(ihi)) then y(ihi)=ytry do 12 j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) 12 continue endif amotry2=ytry return END SUBROUTINE jacobi(a,n,np,d,v,nrot) INTEGER n,np,nrot,NMAX REAL a(np,np),d(np),v(np,np) PARAMETER (NMAX=500) INTEGER i,ip,iq,j REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) do 12 ip=1,n do 11 iq=1,n v(ip,iq)=0. 11 continue v(ip,ip)=1. 12 continue do 13 ip=1,n b(ip)=a(ip,ip) d(ip)=b(ip) z(ip)=0. 13 continue nrot=0 do 24 i=1,50 sm=0. do 15 ip=1,n-1 do 14 iq=ip+1,n sm=sm+abs(a(ip,iq)) 14 continue 15 continue if(sm.eq.0.)return if(i.lt.4)then tresh=0.2*sm/n**2 else tresh=0. endif do 22 ip=1,n-1 do 21 iq=ip+1,n g=100.*abs(a(ip,iq)) if((i.gt.4).and.(abs(d(ip))+ *g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then a(ip,iq)=0. else if(abs(a(ip,iq)).gt.tresh)then h=d(iq)-d(ip) if(abs(h)+g.eq.abs(h))then t=a(ip,iq)/h else theta=0.5*h/a(ip,iq) t=1./(abs(theta)+sqrt(1.+theta**2)) if(theta.lt.0.)t=-t endif c=1./sqrt(1+t**2) s=t*c tau=s/(1.+c) h=t*a(ip,iq) z(ip)=z(ip)-h z(iq)=z(iq)+h d(ip)=d(ip)-h d(iq)=d(iq)+h a(ip,iq)=0. do 16 j=1,ip-1 g=a(j,ip) h=a(j,iq) a(j,ip)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 16 continue do 17 j=ip+1,iq-1 g=a(ip,j) h=a(j,iq) a(ip,j)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 17 continue do 18 j=iq+1,n g=a(ip,j) h=a(iq,j) a(ip,j)=g-s*(h+g*tau) a(iq,j)=h+s*(g-h*tau) 18 continue do 19 j=1,n g=v(j,ip) h=v(j,iq) v(j,ip)=g-s*(h+g*tau) v(j,iq)=h+s*(g-h*tau) 19 continue nrot=nrot+1 endif 21 continue 22 continue do 23 ip=1,n b(ip)=b(ip)+z(ip) d(ip)=b(ip) z(ip)=0. 23 continue 24 continue write(*,*) 'WARNING: too many iterations in jacobi' write(100,*) 'WARNING: too many iterations in jacobi' return END FUNCTION gammq(a,x) REAL a,gammq,x CU USES gcf,gser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.) then write(*,*) 'WARNING: bad arguments in gammq' write(100,*) 'WARNING: bad arguments in gammq' endif if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammq=1.-gamser else call gcf(gammcf,a,x,gln) gammq=gammcf endif return END SUBROUTINE gcf(gammcf,a,x,gln) INTEGER ITMAX REAL a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) CU USES gammln INTEGER i REAL an,b,c,d,del,h,gammln gln=gammln(a) b=x+1.-a c=1./FPMIN d=1./b h=d do 11 i=1,ITMAX an=-i*(i-a) b=b+2. d=an*d+b if(abs(d).lt.FPMIN)d=FPMIN c=b+an/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue write(*,*) 'WARNING: a too large, ITMAX too small in gcf' write(100,*) 'WARNING: a too large, ITMAX too small in gcf' 1 gammcf=exp(-x+a*log(x)-gln)*h return END SUBROUTINE gser(gamser,a,x,gln) INTEGER ITMAX REAL a,gamser,gln,x,EPS PARAMETER (ITMAX=100,EPS=3.e-7) CU USES gammln INTEGER n REAL ap,del,sum,gammln gln=gammln(a) if(x.le.0.)then if(x.lt.0.) then write(*,*) 'WARNING: x < 0 in gser' write(100,*) 'WARNING: x < 0 in gser' endif gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,ITMAX ap=ap+1. del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*EPS)goto 1 11 continue write(*,*) 'WARNING: a too large, ITMAX too small in gser' write(100,*) 'WARNING: a too large, ITMAX too small in gser' 1 gamser=sum*exp(-x+a*log(x)-gln) return END