program PiN c c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c updated 20/6/97 c Main program for Lippman-Schwinger equation desciption c of coupled (piN-piDelta) system. Calculates piN partial wave phase c shifts from threshold to Delta resonance. Also produces the c scattering lengths and differential cross sections. c PiN takes input form factors from chiral quark models c (e.g., chiral color-dielectric model) and searches for best model c parameters (bare masses or glueball mass and coulping constant) c by fitting the nucleon and Delta poles in P11 and P33 channels. c PiN was developed by Dinghui Lu and Rubin Landau in 1994-1995. It c inherits features from the KbarN code (developed by Guangliang c He and Rubin Landau). c Work supported by US DOE, Oregon State University. c*********************************************************************** implicit none c Local variables integer mclock, nch, l, clmax integer ich, jch, ichpt, jchpt, ichpt0, jchpt0 double precision hundred c Shared variables integer Nifty, npts, iopt, n, m, nprint, info, ldfjac integer n_fit, ner, nei, Fidx, Isheet integer Ipvt(3), njev, nfev, mode, maxfev double precision P, W, er1, er2, ei1, ei2 double precision X(3), Fvec(450), Xpole(2) double precision fjac(450, 3), Diag(3), Qtf(3) double precision Wa1(3), Wa2(3), Wa3(3), Wa4(450) double precision factor, ftol, xtol, gtol, epsfcn double precision f0, mgb, alpha, c1, c2, Mbare, mpion, Mass, rcbm double complex zi, Vsuper(128,128,6) double complex Knew(2), enew character*14 reout,imout logical first parameter (zi=(0.0d0, 1.0d0)) parameter (hundred = 1.0d2) common /control/ Nifty(15) c grid.h, gaussian grid points and weights common /grid/ P(64), W(64), npts c bag.h Version 1.4 contains the CCDM parameters. first 3 independent common /bag/ rcbm,f0,mgb,alpha,c1,c2,Mbare(3),Fidx(3) c mass.h Version 1.1 masses of pion, nucleon and delta common /mass/ mpion, Mass(2) c sheet.h: parameters to decide which sheet the energy is on. common /sheet/ Isheet(2) external Chi2, Det1Mvg, mclock data first /.true./ c*********************************************************************** c c start the clock c tic1 = mclock() open(unit=8, file='output',status='unknown') open(unit=22,file='debug', status='unknown') c print out the Banner write(8, 900) write(8, *) 'IBM XL FORTRAN, xlf' c kpinit does the necessary initialization, including open input files call KpInit(er1, er2, ner, ei1, ei2, nei, n_fit, & ftol, xtol, epsfcn) nch = Nifty(2) clmax = Nifty(6) c1111111 if (Nifty(4) .eq. 1) then c calculating the potential only enew = er1+ zi * ei1 write(8, 902) enew do ich = 1, nch Knew(ich) = P(npts) enddo call VBag(npts,P,enew,Knew,Vsuper) c write the potential do 10030 l =1, clmax do jch = 1, nch jchpt0 = (jch-1)*(npts+1) + 1 do ich = jch, nch write(reout, 910) l,ich,jch write(imout, 911) l,ich,jch open(9, file = reout, status='unknown') open(10, file = imout, status='unknown') write(8, 912) l, reout, imout ichpt0 = (ich-1)*(npts+1) + 1 write(9, 916) ((dble(Vsuper(ichpt,jchpt,l)), & ichpt = ichpt0, ichpt0+npts), & jchpt = jchpt0, jchpt0+npts) write(10, 916) ((imag(Vsuper(ichpt,jchpt,l)), & ichpt = ichpt0, ichpt0+npts), & jchpt = jchpt0, jchpt0+npts) close(9) close(10) enddo enddo 10030 continue c222222222 else if (Nifty(4) .eq. 2) then c calculate scattering amplitudes, c phase shifts and scattering lengths once the parameters are fixed if (Nifty(8) .eq. 1) then call Amp(er1, er2, ner, ei1, ei2, nei) else if(Nifty(8) .eq. 2) then call Obs(er1, er2, ner) endif c333333333 else if (Nifty(4) .eq. 3) then c run chi-square fit write(8, *) 'Start to prepare for snls1' iopt = 1 n_fit = 5 m = n_fit + 1 call SetFit(n, X) ldfjac = 450 gtol = 0.0d0 maxfev = 100*(n+1) mode = 1 factor = 1.0d2 nprint = 1 write(*,*)'epsfcn=',epsfcn call snls1(chi2, iopt, m, n, x, Fvec, fjac, ldfjac, ftol, & xtol, gtol, maxfev, epsfcn, Diag, mode, factor, nprint, & info, nfev, njev, Ipvt, Qtf, Wa1, Wa2, Wa3, Wa4) write(8, 920) call PrtPot() write(8, 925) write(8, *) 'INFO', info c Extra run for phashi, PiNNff and p11 c isheet(1) = -1 c isheet(2) = 1 Nifty(4) = 2 Nifty(6) = 6 call Obs(er1,er2,ner) c write(8,*) 'PAREMETERS for POLE search' call PrtPot() Nifty(4) = 4 Nifty(6) = 3 Xpole(1) = 940.-mpion-Mass(1) Xpole(2) = 0. call GetPole(Xpole,xtol,epsfcn) c calc coupling constant er1 = Xpole(1)+mpion+Mass(1) ei1 = Xpole(2) er2 = Mbare(1) ei2 = 0. Nifty(4) = 2 Nifty(6) = 6 Nifty(8) = 1 Nifty(11) = 0 call Amp(er1, er2, ner, ei1, ei2, nei) c444444444 else if (Nifty(4) .eq. 4) then c call snsq searching zero of det(1-VG) for one particular channel c iopt = 2 c n = 2 Xpole(1) = er1-mpion-Mass(1) Xpole(2) = ei1 call GetPole(Xpole,xtol,epsfcn) endif c stop the clock c tic2 = mclock() c on 1 Fortran, one second is 100 tics, not 60 tics as said in manual. c seconds = (tic2 - tic1) * 0.01d0 c write(8, 990) seconds close(8) 900 format(1x, 'PiN --- A Package Calculates PiN-PiDelta system') 902 format('0', 'Energy == ', 2(2x, e10.4)) 910 format('vl+.r',3i1) 911 format('vl+.i',3i1) 912 format('0', 'writing part of Vsuper( , , ', i1, & ') to ', a14, 'and ', a14) 916 format(6(1x,e12.6)) 920 format(/'0',10('='), ' Final Paras + EXTRA OBS and POLE', 10('=')) 925 format(1x, 38('=')) 930 format('0', 10('='), ' Final BS (POLE) Results ', 10('='), /, & 1x, 'Info = ', i4, /, & 1x, 'Re(E) = ', e12.6, 3x, 'Im(E) = ', e12.6, /, & 1x, 'Re(Det) = ', e12.6, 3x, 'Im(Det) = ', e12.6, /, & 1x, 47('=')) 935 format('0', 'Zero not fount. Info = ', i3, /, & 1x, 'See doc of snsq for more details') 990 format('0', 'Total time used: ', f10.3, ' seconds.') end