IAP GITLAB

Commit 0cd8958e authored by Tanguy Pierog's avatar Tanguy Pierog

update branch with EPOS EAS-QGP r7025 using ioclude=4 (effect on air showers)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/branches/EPOS_QGS@7026 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 7ae60eb6
......@@ -49,6 +49,33 @@ fname inihy @CRMCROOT@/tabs/epos.inihy
set pytune 350 !possibility to change PYTHIA tune (for PYTHIA only !)
! modifications for EPOS "QGP"
set nsegce 7
set ptclu 0.
set yradpx 0.5 !0.46
set yradpp 2.5 !0.49 !0.48
!set yradmx 0.12
set yradmi 5.
!set amuseg 40.
set ioclude 4
!set iocova 2
set amdrmin 7.
set amdrmax 70.
set reminv 0.15
set rstras(2) 0.2
set wgtdiq 0.11 !22
set wgtqqq(1) 0.1 !22
set wgtqqq(2) 0.1 !22
set wgtqqq(3) 0.1 !22
set alpdi(1) 0.95
!set alpdi(1) 0.8
!set alpndi(1) 1.5
set rexres(1) 0.!15
!!ImpactParameter
!set bminim 0 !works with epos
!set bmaxim 4
......
......@@ -426,7 +426,8 @@ c------------------------------------------------------------------------
common/hadr8/alpqua,alppar,alpsea,alpval,alpdiq,alplea(4),alpdif
real alpndi,alpdi,ptsendi,zdrinc,zmsinc,ptsems
integer irzptn
common/hadr14/alpndi,alpdi,ptsendi,zdrinc,zmsinc,ptsems,irzptn
common/hadr14/alpndi(2),alpdi(2),ptsendi,zdrinc,zmsinc,ptsems
& ,irzptn
real zbcut,zopinc,zipinc,zoeinc,xmxrem
common/hadr15/zbcut,zopinc,zipinc,zoeinc,xmxrem
real fkainc,fkamax,zodinc,zbrmax,zdfinc,xzcut,ptvpom
......
......@@ -10,7 +10,7 @@ c-----------------------------------------------------------------------
include 'epos.incsem'
common/record/maxrec(2),irecty(30,2)
common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi/cifset/ifset
common /ems12/iodiba,bidiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
c &,fndpmjet,fndpmjetpho,fndpmpath !qgs-II????????
c common/dpmjetfname/ fndpmjet,fndpmjetpho,fndpmpath
......@@ -393,8 +393,10 @@ c ... up to here.
rexres(3)=0.5 !kaon remnant excitation probability in nucleus
rexres(4)=1. !charm remnant excitation probability in nucleus
alpdif=0.45 !alpha mass diffractive for cross section and metropolis
alpdi=0.45 !alpha mass diffractive
alpndi=1.6 !alpha mass nondiffractive
alpdi(1)=0.45 !alpha mass diffractive
alpdi(2)=0.45 !alpha mass diffractive
alpndi(1)=1.6 !alpha mass nondiffractive
alpndi(2)=1.6 !alpha mass nondiffractive
alpsea=0.3 !alpha string end x for sea parton
alpval=0.3 !alpha string end x for valence parton
alpdiq=0.3 !alpha string end x for sea diquark
......@@ -407,8 +409,8 @@ c ... up to here.
alpdro(2)=0.3 !pt of leading droplet
alpdro(3)=1.6 !alpha mass of leading droplet (iept=3)
iodiba=0. ! if iodiba=1, study H-Dibaryon (not used (see ProRef in epos-ems ????)
bidiba=0.030 ! epsilon of H-Dibaryon
iodiba=0. ! if iodiba=1, study H-Dibaryon
bidiba=0.20 ! epsilon of H-Dibaryon
c screening splitting +++
......@@ -457,7 +459,7 @@ c masses
amhadr(7)=.548 !eta mass
amhadr(8)=3.5 !J/psi mass
c qcmass=1.6 !c-quark mass (in idmass = 1.2 )
c amhdibar=2.200 !h-dibaryon mass !not used any more
amhdibar=1.80 !h-dibaryon mass !not used any more
qmass(0)=pmqq !diquark effective bounding energy (for pt distribtions)
isospin(0)=0
call idmass(1,qumass)
......@@ -932,6 +934,7 @@ c-----------------------------------------------------------------------
infragm=2 !nuclear fragmentation
iospec=18 !option for particle species (allow deuterium by default)
c cross section based parameters
......@@ -976,8 +979,10 @@ c remnant excitation
rexres(2)=0. !nucleon remnant excitation probability in nucleus
rexres(3)=0.15 !kaon remnant excitation probability in nucleus
alpdif=0.7 !alpha mass diffractive for cross section and metropolis
alpdi=1.05 !alpha mass diffractive
alpndi=2. !1.65 !alpha mass nondiffractive
alpdi(1)=1.05 !alpha mass diffractive
alpdi(2)=1.05 !alpha mass diffractive
alpndi(1)=2. !1.65 !alpha mass nondiffractive
alpndi(2)=2. !1.65 !alpha mass nondiffractive
alpdro(3)=2.5 !alpha mass of leading droplet (iept=3)
alpdro(2)=1.5 !alpha mass of leading droplet (iept=3)
zmsinc= 0. !increase of remant minimum mass and decrease alpha (increase remnant mass with iept=3)
......@@ -2996,6 +3001,7 @@ c-----------------------------------------------------------------------
double precision rcproj,rctarg
common/geom1/rcproj,rctarg
common/photrans/phoele(4),ebeam
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
external sptj
......@@ -3493,7 +3499,7 @@ c---------------------------------------------------------------------
data nappl /0/
common/record/maxrec(2),irecty(30,2)
common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi
common /ems12/iodiba,bidiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
c &,fndpmjet,fndpmjetpho,fndpmpath
c common/dpmjetfname/ fndpmjet,fndpmjetpho,fndpmpath
......@@ -4420,7 +4426,7 @@ c unified parameters
if(linex(ix:jx).eq.'sloreg')sloreg=sngl(val)
if(linex(ix:jx).eq.'gamreg')gamreg=sngl(val)
if(linex(ix:jx).eq.'r2reg' )r2reg= sngl(val)
c if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
if(linex(ix:jx).eq.'ptdiff')ptdiff=sngl(val)
if(linex(ix:jx).eq.'alppar')alppar=sngl(val)
if(linex(ix:jx).eq.'alpsea')alpsea=sngl(val)
......@@ -4431,8 +4437,10 @@ c if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
if(linex(ix:jx).eq.'alplea(3)')alplea(3)=sngl(val)
if(linex(ix:jx).eq.'alplea(4)')alplea(4)=sngl(val)
if(linex(ix:jx).eq.'alpdif')alpdif=sngl(val)
if(linex(ix:jx).eq.'alpdi')alpdi=sngl(val)
if(linex(ix:jx).eq.'alpndi')alpndi=sngl(val)
if(linex(ix:jx).eq.'alpdi(1)')alpdi(1)=sngl(val)
if(linex(ix:jx).eq.'alpdi(2)')alpdi(2)=sngl(val)
if(linex(ix:jx).eq.'alpndi(1)')alpndi(1)=sngl(val)
if(linex(ix:jx).eq.'alpndi(2)')alpndi(2)=sngl(val)
if(linex(ix:jx).eq.'ammsqq')ammsqq=sngl(val)
if(linex(ix:jx).eq.'ammsqd')ammsqd=sngl(val)
if(linex(ix:jx).eq.'ammsdd')ammsdd=sngl(val)
......@@ -5613,7 +5621,7 @@ c elastic event
if(iret.gt.0)goto 3
maxfra=nptl
if(ish.ge.2.and.model.eq.1)
& call alist('list after fragmentation&',nptlx,nptl)
& call alist('list after fragmentation&',1,nptl)
if(irescl.eq.1)then
call utghost(iret)
if(iret.gt.0)goto 3
......
......@@ -772,7 +772,7 @@ c writes /cptl/
c-----------------------------------------------------------------------
include "epos.inc"
include "epos.incems"
double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
common/nucl3/phi,bimp
common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
......@@ -1316,132 +1316,6 @@ c-----------------------------------------------------------------------
end
c----------------------------------------------------------------------
subroutine geoglauber
c----------------------------------------------------------------------
include 'epos.inc'
include 'epos.incems'
common/geom/rmproj,rmtarg,bmax,bkmx
common/nucl3/phi,bimp
common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
* ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
dimension kopj(mamx),kotg(mamx)
logical phi_zero
phi_zero=.false.
! calculate collision number according to Glauber ------------------
c rs=r2had(iclpro)+r2had(icltar)+slopom*log(sy)
c bglaub2=4.*.0389*rs
xs=sigine
bglaub=max(1.,xs/10./pi) !10= fm^2 -> mb
bglaub=sqrt(bglaub)
nglevt=0
do i=1,maproj
kopj(i)=0
enddo
do i=1,matarg
kotg(i)=0
enddo
if(abs(phimin).lt.1e-5.and.abs(phimax).lt.1e-5)phi_zero=.true.
if(phi_zero)then
x2av=0
y2av=0
xyav=0
xav=0
yav=0
bx=cos(phi)*bimp
by=sin(phi)*bimp
else
bx=0
by=0
endif
do ko=1,koll
ip=iproj(ko)
it=itarg(ko)
r=bk(ko)
if(phi_zero)then
r2=(xproj(ip)+bx-xtarg(it))**2+(yproj(ip)+by-ytarg(it))**2
if((maproj.gt.1.or.matarg.gt.1).and.abs(r**2-r2).gt.1e-3)then
print*,'geoglauber: r**2,r2:',r**2,r2
stop
endif
endif
if(r.le.bglaub)then
nglevt=nglevt+1
kopj(ip)=kopj(ip)+1
kotg(it)=kotg(it)+1
endif
enddo
ng1evt=0
ng2evt=0
do i=1,maproj
if(kopj(i).ge.1)ng1evt=ng1evt+1
if(kopj(i).ge.2)ng2evt=ng2evt+1
enddo
do i=1,matarg
if(kotg(i).ge.1)ng1evt=ng1evt+1
if(kotg(i).ge.2)ng2evt=ng2evt+1
enddo
if(phi_zero)then
do i=1,maproj
if(kopj(i).ge.1)then
x=xproj(i)+bx/2
y=yproj(i)+by/2
xav=xav+x
yav=yav+y
x2av=x2av+x*x
y2av=y2av+y*y
xyav=xyav+x*y
!print*,x,y
endif
enddo
do i=1,matarg
if(kotg(i).ge.1)then
x=xtarg(i)-bx/2
y=ytarg(i)-by/2
xav=xav+x
yav=yav+y
x2av=x2av+x*x
y2av=y2av+y*y
xyav=xyav+x*y
!print*,x,y
endif
enddo
if(ng1evt.gt.0)then
x2av=x2av/float(ng1evt)
y2av=y2av/float(ng1evt)
xyav=xyav/float(ng1evt)
xav= xav /float(ng1evt)
yav= yav /float(ng1evt)
eglevt=(y2av-x2av)/(y2av+x2av)
fglevt=sqrt( (y2av-x2av)**2+4*(xyav-xav*yav)**2 ) /(y2av+x2av)
rglevt=ng2evt/float(ng1evt)
sglevt=pi*sqrt(x2av)*sqrt(y2av)
!~~~~~~~~~~~~~~~~~~~~
xx=x2av ! <x**2>
yy=y2av ! <y**2>
xy=xyav ! <x*y>
dta=0.5*abs(xx-yy)
ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
yy=ev1
xx=ev2
fglevt=(yy-xx)/(yy+xx)
!~~~~~~~~~~~~~~~~~~~~
endif
!if(ng1evt.gt.0)print*,bimp,eglevt,bx,by,phi
else
eglevt=-1.
fglevt=-1.
rglevt=-1.
sglevt=-1.
endif
end
......
......@@ -128,7 +128,10 @@ c-----------------------------------------------------------------------
common/cspez4/ffstat(2,0:mspez+2) /ctfo/tfo
real u(3)
io3=ioclude-3
if(io3.lt.1.or.io3.gt.2)stop'in GraCan: wrong ioclude (140808) '
if(io3.lt.1.or.io3.gt.2)then
print *,'ioclude', ioclude
stop'in GraCan: wrong ioclude (140808) '
endif
yie= volu * ffstat(io3,nspez) * ffstat(io3,mspez+2)
np=yie
if(rangen().le.yie-np)np=np+1
......@@ -1490,7 +1493,7 @@ c--------------------------------------------------------------------
!-----------------------------------------------------------
include 'epos.inc'
include 'epos.inchy'
parameter (mspez=54,msto=2)
parameter (mspez=54,msto=3)
common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
common/cspez2/kspez,jspez(mspez)
common/cspez3/fuga(mspez)
......@@ -1499,7 +1502,9 @@ c--------------------------------------------------------------------
common/sp1/ffstatSave(2,0:mspez+2,msto),fugaSave(mspez,msto)
common/sp2/iocludeSave(msto),kspezSave(msto),nspezSave(msto)
common/sp3/tfoSave(msto),meosmuSave(msto),epscri3Save(msto)
if(isens.eq.1)then !~~~~~~~~~~~~~~store~~~~~~~~~~~~~~
c write(ifch,*)'mani',isens,isto,ioclude,iocludeSave(isto)
if(isens.eq.1)then !~~~~~~~~~~~~~~store~~~~~~~~~~~~~~
iocludeSave(isto)=ioclude
kspezSave(isto)=kspez
nspezSave(isto)=nspez
......@@ -1542,23 +1547,24 @@ c--------------------------------------------------------------------
ioclude=5 !choice of hadronization (4 or 5)
kspez=1 !choice of
nspez=54 ! particles
tfo=0.130 !freeze out temperature used in this routine
tfo=0.1402 !0.1176 !0.130 !freeze out temperature used in this routine
meosmu=0 !no chemical potentials used in this routine
epscri(3)=0.0765 !should be consitent with tfo see EoS table below
epscri(3)=0.125 !0.04 !0.0765 !should be consitent with tfo see EoS table below
call DefineParticles
elseif(isens.eq.93)then !~~~~~~~~~~~~~~redefine set3~~~~~~~~~~~~~
!used for application micro
ioclude=5 !choice of hadronization (4 or 5)
kspez=1 !choice of
nspez=54 ! particles
tfo=0.130 !freeze out temperature used in this routine
tfo=0.12157 !0.130 !freeze out temperature used in this routine
meosmu=0 !no chemical potentials used in this routine
epscri(3)=0.0765 !should be consitent with tfo see EoS table below
epscri(3)=0.05 !0.0765 !should be consitent with tfo see EoS table below
call DefineParticles
endif
!--------------------------!
! epscri(3) ! tfo !
!--------------------------!
! 0.045 ! 0.11954 !
! 0.045 ! 0.11954 ! less strange baryons
! 0.050 ! 0.12157 !
! 0.055 ! 0.12327 !
! 0.060 ! 0.12496 !
......@@ -1574,7 +1580,7 @@ c--------------------------------------------------------------------
! 0.110 ! 0.13749 !
! 0.115 ! 0.13851 !
! 0.120 ! 0.13952 !
! 0.125 ! 0.1402 !
! 0.125 ! 0.1402 ! more strange baryons
!--------------------------!
end
......
This diff is collapsed.
This diff is collapsed.
......@@ -25,7 +25,7 @@ c computes charge of particle with ident code id
c ichrg must be dimensioned nqlep+12
c-----------------------------------------------------------------------
dimension ichrg(53),ifl(3)
data ichrg/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,1,1,2,2*0
data ichrg/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,3,3,3,2*0
*,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,0,-3,3,0
*,3,0,0,0,3,3,3,6,6,6,0/
idabs=iabs(id)
......@@ -315,7 +315,7 @@ c quarks, leptons, etc
ifl3=0
jspin=0
index=idabs
if(idabs.lt.20) return
if(idabs.lt.20.or.idabs.eq.30) return
c define index=20 for ks, index=21 for kl
index=idabs+1
if(id.eq.20) index=20
......@@ -546,6 +546,7 @@ c-----------------------------------------------------------------------
c returns the mass of the particle with ident code id.
c (deuteron, triton and alpha mass come from Gheisha ???)
c-----------------------------------------------------------------------
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
dimension ammes0(15),ammes1(15),ambar0(30),ambar1(30)
dimension amlep(52)
parameter ( nqlep=41,nmes=2)
......@@ -600,10 +601,10 @@ c 3/2+ baryon mass table
c entry
id=idi
amass=0.
ctp060829 if(iabs(id).eq.30)then
ctp060829 amass=amhdibar
ctp060829 return
ctp060829 endif
if(iabs(id).eq.30)then
amass=amhdibar
return
endif
if(idi.eq.0)then
id=1120 !for air target
elseif(abs(idi).ge.1000000000)then
......@@ -1741,6 +1742,10 @@ c-----------------------------------------------------------------------
ic(1)=222000
ic(2)=000000
return
elseif(id.eq.-30)then
ic(1)=000000
ic(2)=222000
return
endif
if(iabs(id).lt.1e8)then
ix=iabs(id)/10
......
This diff is collapsed.
......@@ -104,12 +104,11 @@ c----------------------------------------------------------------------
! calculate collision number according to Glauber -----------------------
bglaub=sqrt(sigine/10./pi)
c nglevt=0
c do ko=1,koll
c r=bk(ko)
c if(r.le.bglaub)nglevt=nglevt+1
c enddo
call geoglauber
nglevt=0
do ko=1,koll
r=bk(ko)
if(r.le.bglaub)nglevt=nglevt+1
enddo
! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)-----------
......
......@@ -171,7 +171,9 @@ c -----
do i=nptlpt+1,nptl
if(mod(istptl(i),10).eq.0)then
pptl(1,i)=pptl(1,i)-sngl(pt1soll)
c . +2.*max(0.,pptl(1,i)-0.5)
pptl(2,i)=pptl(2,i)-sngl(pt2soll)
c . +2.*max(0.,pptl(2,i)-0.5)
pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2
& +pptl(3,i)**2+pptl(5,i)**2)
endif
......@@ -335,7 +337,7 @@ c mark ghosts
c -----------
do j=nptlpt+1,nptl
if(istptl(j).le.1.and.pptl(4,j).gt.0.d0)then
if(istptl(j).eq.0)then !Don't fix mass of decayed particles (to keep width)
if(iLHC.eq.0.or.mod(abs(idptl(j)),10).le.1)then !for LHC tune don't fix mass of resonnances (to keep width)
amass=pptl(5,j)
call idmass(idptl(j),amass)
if(abs(idptl(j)).gt.100.and.
......@@ -1561,8 +1563,6 @@ c-----------------------------------------------------------------------
parameter (nflav=6)
integer jcp(nflav,2),jcm(nflav,2)
imod=5 !4 for resonance and 5 for minimum mass
do i=1,nflav
do j=1,2
if(jcp(i,j).ne.0)goto1
......@@ -1574,7 +1574,7 @@ c-----------------------------------------------------------------------
kec=jcm(4,1)-jcm(4,2)
keb=jcm(5,1)-jcm(5,2)
ket=jcm(6,1)-jcm(6,2)
utamnx=utamnu(keu,ked,kes,kec,keb,ket,imod)
utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
return
1 continue
......@@ -1589,7 +1589,7 @@ c-----------------------------------------------------------------------
kec=jcp(4,1)-jcp(4,2)
keb=jcp(5,1)-jcp(5,2)
ket=jcp(6,1)-jcp(6,2)
utamnx=utamnu(keu,ked,kes,kec,keb,ket,imod)
utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
return
2 continue
......@@ -1602,16 +1602,16 @@ c-----------------------------------------------------------------------
ke=keu+ked+kes+kec+keb+ket
if(mod(ke+1,3).eq.0)then
keu=keu+1
amms1=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
keu=keu-1
ked=ked+1
amms2=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
elseif(mod(ke-1,3).eq.0)then
keu=keu-1
amms1=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
keu=keu+1
ked=ked-1
amms2=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
else
amms1=0
amms2=0
......@@ -1628,16 +1628,16 @@ c-----------------------------------------------------------------------
ke=keu+ked+kes+kec+keb+ket
if(mod(ke+1,3).eq.0)then
keu=keu+1
amms3=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
keu=keu-1
ked=ked+1
amms4=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
elseif(mod(ke-1,3).eq.0)then
keu=keu-1
amms3=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
keu=keu+1
ked=ked-1
amms4=utamnu(keu,ked,kes,kec,keb,ket,imod)
amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
else
call utstop('utamnx: no singlet possible (2)&')
endif
......@@ -3675,6 +3675,9 @@ c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
function fremnux(jc)
c-----------------------------------------------------------------------
real bidiba,amhdibar
integer iodiba
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
real pnll,ptq
common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
parameter (nflav=6)
......@@ -3682,6 +3685,15 @@ c-----------------------------------------------------------------------
c ic(1)=210000
c ic(2)=0
c call iddeco(ic,jc)
if(iodiba.eq.1)then
if((jc(1,1).eq.2.and.jc(2,1).eq.2.and.jc(3,1).eq.2
. .and.jc(1,2)+jc(2,2)+jc(3,2).eq.0)
. .or.(jc(1,2).eq.2.and.jc(2,2).eq.2.and.jc(3,2).eq.2
. .and.jc(1,1)+jc(2,1)+jc(3,1).eq.0))then
fremnux=amhdibar-bidiba
return
endif
endif
keu=jc(1,1)-jc(1,2)
ked=jc(2,1)-jc(2,2)
kes=jc(3,1)-jc(3,2)
......
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment