Текст программы modmd82

program dat(input,output);

label ll; const ks=12; ps=3; ProgID=82; r11=0.600;r12=0.240;s=3;

var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1,su1,sum1,su2,sum2:real;

x2,y2,z2,vx00,vyz00,vz00,vy00,Tm,Tmm,r0,r1,r2,ugol,ac,ugol_x,ugol_y,Tsr,Tss,Vsr,VVsr:real;

vxTmm,vxTmmv,vyzTmm,vyTmm,vyTmmv,vzTmm,vzTmmv,tid,Tsss,NNN, R,T,L,rA,r01,r02,kpr:real;

j,p,m1,k,kR,kL,f3,f4,sk:integer; n,i,n1,n2:longint; isp_No,m3,tekp:integer; defl:boolean;

inp1,out2,out3,out4:text; Fdir,conc,Nup,Ndn:real;

nk: array[1..2*ks,1..ps] of real;

nk1: array[1..2*ks,1..ps] of real; nc: array[1..ps,1..2] of real; d: array [1..ps] of real;

procedure Maxwell;

begin

repeat

v:=(random(5)+random);

y:=v*v*exp(1-v*v);

y0:=random;

until y0<=y;

if (v<=0.05) then v:=0.05;

end;

procedure scattering(Mx,My:real);

var vx0,vyz0,vy0,vz0,gamma,beta,alfa:real;

begin

repeat

beta:=pi/2*random;

y:=2*sin(beta)*cos(beta);

y0:=random;

until y0<=y;

Maxwell;

vv:=(1-ac)*vv*vv+ac*v*v+ac*(1-ac)*vv*vv*((v*v)/1.5-1);

if (vv<=0.0025) then vv:=0.05

else vv:=sqrt(vv);

alfa:=2*pi*random;

vx0:=vv*cos(beta); vyz0:=vv*sin(beta);

vy0:=vyz0*cos(alfa); vz0:=vyz0*sin(alfa);

if Mx<2 then begin vx:=Mx*vx0; vy:=vy0; vz:=Mx*vz0;end

else if not(y2*z2=0) then begin

if (y2*My)>0 then gamma:=arctan(z2/y2)+Pi-arctan(vy0/vx0)

else gamma:=arctan(z2/y2)-arctan(vy0/vx0);

vx:=-vz0; vy:=sqrt(vx0*vx0+vy0*vy0)*cos(gamma);

vz:=sqrt(vx0*vx0+vy0*vy0)*sin(gamma); end

else if y2=0 then begin

vx:=-vz0;vy:=-vy0*My*z2/abs(z2);vz:=-vx0*My*z2/abs(z2);end

else begin vx:=-vz0;vy:=-vx0*My*y2/abs(y2);vz:=vy0*My*y2/abs(y2);end;

end;

procedure counttrack;

label 2,3;

var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod:real;

f0,f1,f2,fr1,fr2,fa2:real;

procedure countmolec;

var w: 0..ps; b1,b2,b3,b4,a1,a2:boolean;

d1,d2,d3,d4,tau1,tau2,tau3,tau4,tkon1,tkon2,t0,t00:real; label ii;

begin

Tsr:=Tsr+abs(tmax-tmin);

for w:=0 to ps-1 do

begin

tkon1:=((L*w/ps)-x1)/vx; tkon2:=((L*(w+1)/ps)-x1)/vx;

p:=w+1;

a1:=((tkon1>=tmin)and(tkon1<=tmax));

a2:=((tkon2>=tmin)and(tkon2<=tmax));

if a1 then if a2 then if vx>0 then begin t0:=tkon1;t00:=tkon2;end

else begin t0:=tkon2; t00:=tkon1;end

else if vx>0 then begin t0:=tkon1; t00:=tmax;end

else begin t0:=tmin; t00:=tkon1;end

else if a2 then if vx>0 then begin t0:=tmin; t00:=tkon2;end

else begin t0:=tkon2; t00:=tmax;end

else if (tkon1<tmin)and(tkon2>tmax) or (tkon1>tmax)and(tkon2<tmin)

then begin t0:=tmin; t00:=tmax;end

else goto ii;

d1:=y1-z1*vy/vz; d2:=d1;

tau1:=(d1-y1)/vy; tau2:=(d2-y1)/vy;

for kR:=1 to ks do

begin

d3:=vz*y1-vy*z1;

d4:=d3/(vz*cos(kR*Pi/ks)+vy*sin(kR*Pi/ks));

d3:=d3/(vz*cos(kR*Pi/ks)-vy*sin(kR*Pi/ks));

tau3:=(d3*cos(kR*Pi/ks)-y1)/vy; tau4:=(d4*cos(kR*Pi/ks)-y1)/vy;

b1:=((d1>0)and(tau1<=t00)and(tau1>=t0));

b2:=((d2>0)and(tau2<=t00)and(tau2>=t0));

b3:=((d3>0)and(tau3<=t00)and(tau3>=t0));

b4:=((d4>0)and(tau4<=t00)and(tau4>=t0));

if b1 then if b3 then nk[kR,p]:=nk[kR,p]+abs(tau1-tau3)

else if ((tau3-tau1)*d3)>0 then nk[kR,p]:=nk[kR,p]+abs(tau1-t00)

else nk[kR,p]:=nk[kR,p]+abs(tau1-t0)

else if b3 then if ((tau1-tau3)*d1)>0 then nk[kR,p]:=nk[kR,p]+abs(tau3-t00)

else nk[kR,p]:=nk[kR,p]+abs(tau3-t0)

else if (((d1>0)and(d3>0)and(((tau1>t00)and(tau3<t0))or((tau1<t0)and(tau3>t00))))

or ((d1*d3<0)and(((d3*(tau1-tau3)>0)and(tau1>t00)and(tau3>t00))

or((d1*(tau3-tau1)>0)and(tau1<t0)and(tau3<t0)))))

then nk[kR,p]:=nk[kR,p]+abs(t00-t0);

kL:=(2*ks+1)-kR;

if b2 then if b4 then nk[kL,p]:=nk[kL,p]+abs(tau2-tau4)

else if ((tau4-tau2)*d4)>0 then nk[kL,p]:=nk[kL,p]+abs(tau2-t00)

else nk[kL,p]:=nk[kL,p]+abs(tau2-t0)

else if b4 then if ((tau2-tau4)*d2)>0 then nk[kL,p]:=nk[kL,p]+abs(tau4-t00)

else nk[kL,p]:=nk[kL,p]+abs(tau4-t0)

else if (((d2>0)and(d4>0)and(((tau2>t00)and(tau4<t0))or((tau2<t0)and(tau4>t00))))

or ((d2*d4<0)and(((d4*(tau2-tau4)>0)and(tau2>t00)and(tau4>t00))

or((d2*(tau4-tau2)>0)and(tau2<t0)and(tau4<t0)))))

then nk[kL,p]:=nk[kL,p]+abs(t00-t0);

tau1:=tau3;tau2:=tau4;d1:=d3;d2:=d4;

end;

ii:end;

end;

begin

tmax:=0; tmin:=0;

x1:=x2;y1:=y2;z1:=z2;

f0:=vy*vy+vz*vz; f1:=(vy*y1+vz*z1)/f0;

f2:=f0-sqr(vz*y1-vy*z1);

fa2:=rA*rA*f0-sqr(vz*y1-vy*z1);

fr1:=r1*r1*f0-sqr(vz*y1-vy*z1);

fr2:=r2*r2*f0-sqr(vz*y1-vy*z1);

if f2>0 then

begin

f2:=sqrt(f2)/f0; taumin:=-f1-f2; taumax:=-f1+f2;

if not(vx=0) then if vx>0 then begin taukon:=(L-x1)/vx; f3:=-1;end

else begin taukon:=(0-x1)/vx; f3:=1;end

else taukon:=(2/vyz)+0.1;

tAnod:=-1;

if fr2>0 then

begin

fr2:=sqrt(fr2)/f0; tmin2:=-f1-fr2; tmax2:=-f1+fr2;

if fr1>0 then

begin

fr1:=sqrt(fr1)/f0; tmin1:=-f1-fr1; tmax1:=-f1+fr1;

if fa2>0 then begin fa2:=sqrt(fa2)/f0; tAnod:=-f1-fa2; end;

if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2

else if (tmin2<=0)and(tmin1>0) then Tmin:=0

else goto 2;

if (tmin1<taukon) then Tmax:=tmin1

else Tmax:=taukon;

countmolec;

2: if (tmax1>0)and(tmax1<taukon) then Tmin:=tmax1

else if (tmax1<=0)and(tmax2>0) then Tmin:=0

else goto 3;

if (tmax2<taukon) then Tmax:=tmax2

else Tmax:=taukon;

if (tAnod<0) then countmolec;

goto 3;

end;

if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2

else if (tmin2<=0)and(tmax2>0) then Tmin:=0

else goto 3;

if (tmax2>taukon) then Tmax:=taukon

else Tmax:=tmax2;

countmolec;

end;

3:end;

if taumax<taukon then begin taukon:=taumax; f3:=2;f4:=1;end;

if (tAnod>0)and(tAnod<taukon) then begin taukon:=tAnod; f3:=2;f4:=-1;end;

x2:=x1+taukon*vx; y2:=y1+taukon*vy; z2:=z1+taukon*vz;

Tss:=Tss+taukon/tid;Tsss:=Tsss+taukon; {ў Tss ­ Є Ї«Ёў Ґ taukon/tid; Tsss - бв ஥ §­ 祭ЁҐ Tss}

if x2>L then x2:=L; if x2<0 then x2:=0;

end;

procedure byben;

var y0,y,sum:real;

begin

m3:=0;

if(r0>r11-r12)and(r0<r11+r12) then

begin

sum:=(r11+r12+r0)/2;

y0:=1/3+4/pi*2*arctan(sqrt((sum-r11)*(sum-r12)*(sum-r0)/sum)/(sum-r12));

y:=random;

if(y<=y0) then m3:=1;

end;

end;

{procedure vvod (var Nsu:mas);

begin

assign(out4,'c:\80Nsum.txt');

reset(out4);

for i1:=1 to s do

for j1:=1 to m do

read(out4,Nsu[i1,j1]);

{writeln(a[1,2]:7:4);

writeln(a[2,2]:7:4);

writeln(a[3,2]:7:4);

close(out4);

end;}

function summa:real;

var su1,sum1,su2,sum2,Koef:real;

begin

su1:=0;

for p:=1 to s do

su1:=su1+nc[p,1];

sum1:=su1/3;

su2:=0;

for p:=1 to s do

su2:=su2+nc[p,2];

sum2:=su2/3;

writeln(out2,'Nup=':8,'Ndn=':8);

writeln(out2,sum1:8:4,sum2:8:4);

Koef:=(sum2-sum1)/(sum1+sum2);

writeln(out2,'Koef=':8);

writeln(out2,Koef:8:4);

end;

begin

Randomize;

Assign(inp1,'c:\80inp.txt'); Reset(inp1);

writeln; writeln('Start, program ID=',ProgID:2); isp_No:=0;

repeat

Assign(out2,'c:\80all.txt'); Append(out2);

isp_No:=isp_No+1;

Read(inp1,n,L,rA,r01,r02,ugol,kpr,r1,r2,v0,ugol_x,ugol_y,R,Tm,Tmm,T,ac,m1);

writeln('---------------------------€б室­лҐ-¤ ­­лҐ------------------ь-ў-бҐаЁЁ:',isp_No:2);

writeln('n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);

writeln(n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);

writeln('Vo=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);

writeln(v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);

w0:=v0/sqrt(R*T); vTmm:=sqrt(Tmm/T); Fdir:=n;

vTm:=sqrt(Tm/T);vn:=0;Tsr:=0;Tss:=0;Vsr:=0;VVsr:=0;n1:=0;n2:=0;NNN:=0;

vx00:=w0*cos(ugol_x/180*pi);vyz00:=w0*sin(ugol_x/180*pi);

vy00:=vyz00*cos(ugol_y/180*pi);vz00:=vyz00*sin(ugol_y/180*pi);

vxTmm:=vTmm*cos(ugol_x/180*pi);vyzTmm:=vTmm*sin(ugol_x/180*pi);

vyTmm:=vyzTmm*cos(ugol_y/180*pi);vzTmm:=vyzTmm*sin(ugol_y/180*pi);

a:=(4*r12*r12 +(4/3)*r11*r12)*ugol/360;

q:=(r2*r2-r1*r1)*L/ps/(2*ks)/a;

q1:=(r2*r2-r1*r1)/ps/(2*ks)/a*(1-rA*rA); tekp:=0; {q1=(q/L)*(1-rA*rA), (1-rA^2) - Ї«®й ¤м Ґ¦н«ҐЄвத­®© з бвЁ в®а楢}

for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=0;

for i:=1 to n do

begin

if ((i/n*100) > (tekp+1)) then begin write(chr(8),chr(8),chr(8),chr(8),'%',tekp:3); tekp:=tekp+1 end;

Maxwell; defl:=false;

vxTmmv:=v*vxTmm;vyTmmv:=v*vyTmm;vzTmmv:=v*vzTmm;

Maxwell; v:=v*vTm+0.0000000001;

{teta:=pi*random;} teta:=2*random-1;

{vx:=v*cos(teta)+vx00+vxTmmv;} vx:=v*teta+vx00+vxTmmv;

if not(vx=0) then begin

if(vx>0) then begin sk:=-1;x2:=0; n1:=n1+1 end

else begin sk:=1; x2:=L; n2:=n2+1 end;

repeat y2:=2*Random-1; z2:=2*Random-1;

r0:=sqrt(y2*y2+z2*z2);

byben;

{xkpr:=Random;}

until (y2*y2+z2*z2<r02*r02)and(y2*y2+z2*z2>r01*r01)and(m3=1)and(((sk=1)and(-cos(pi*ugol/360)+y2/sqrt(y2*y2+z2*z2)>0))

or((sk=-1)and(cos(pi-pi*ugol/360)-y2/sqrt(y2*y2+z2*z2)>0)));

{vyz:=v*sin(teta);} vyz:=v*sqrt(1-teta*teta);

ksi:=2*pi*random;

vy:=vyz*cos(ksi)+vy00+vyTmmv; vz:=vyz*sin(ksi)+vz00+vzTmmv;

vv:=sqrt(vx*vx+vy*vy+vz*vz); VVsr:=VVsr+vv;

vn:=vn+abs(vx); j:=0; tid:=abs(L/vx); {<- (!) ¤«п Ї®«­. ®вЄа. вагЎл, ®ЇҐа в®а ЇаЁбў. ¤«п tid в®«мЄ® §¤Ґбм!}

for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=0;Tsss:=0;

repeat counttrack;

r0:=sqrt(y2*y2+z2*z2);

xkpr:=Random;

if (xkpr<kpr)and(abs(f3)=1)and(y2*y2+z2*z2>r01*r01)and(y2*y2+z2*z2<r02*r02)and(m3=1)and(((f3=-1)

and(-cos(pi*ugol/360)+y2/sqrt(y2*y2+z2*z2)>0))

or((f3=1)and(cos(pi-pi*ugol/360)-y2/sqrt(y2*y2+z2*z2)>0)))

then goto ll {Ґб«Ё Ї®Ї «ў §®­г Їа®§а з­®бвЁ, в® ўл«Ґв}

else begin scattering(f3,f4); defl:=true end {Ё­ зҐ, ®ва ¦Ґ­ЁҐ}

until j=1;

ll: Vsr:=Vsr+vv; if (defl) then Fdir:=Fdir-1;

for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=nk[k,p]/tid; { ®в­®иҐ­ЁҐ ўаҐҐ­ ॠ«.(¤«п Њ€Џ) Ё Ё¤Ґ «. (¤«п ®вЄалв®© вагЎл)}

for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=nk1[k,p]+nk[k,p];

end;

end;

for k:=1 to 2*ks do for p:=1 to ps do NNN:=NNN+nk1[k,p]; { cгЁа㥠ўбҐ ®в­®иҐ­Ёп ўаҐҐ­ ॠ«. Ё Ё¤Ґ «. Ї® з.®ЎкҐ }

ww0:=vn/n; Tss:=Tss*a/n*(1-rA*rA);NNN:=NNN/n/q1/ps/(2*ks); conc:=Tsr/n*ww0/q/ps/(2*ks); Fdir:=Fdir/n;

{conc - бв ஥ §­ з. ®в­.Є®­ж.ў Њ€Џ; ў Tss ¤® бЁе Ї®а ­ Є Ї«Ёў «®бм ®в­®иҐ­ЁҐ taukon/tid, ⥯Ґам Tss ®в­®аЁа.; ­®аЁа㥠NNN}

write(chr(8),chr(8),chr(8),chr(8));

writeln('------------------------------ђҐ§г«мв вл--------------------------------');

writeln('N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vўе|=':8,'|Vўле|=':8,'Tss=':8,'Fdir=':7);

writeln(conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:7:3);

writeln('------------------------------------------------------------------------');

writeln(out2,'ProgID=':8,'n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);

writeln(out2,ProgID:8,n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);

writeln(out2,'V0=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);

writeln(out2,v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);

writeln(out2,'N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vin|=':8,'|Vout|=':8,'Tss=':8,'Fdir=':8);

writeln(out2,conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:8:3);

writeln(out2,'Nup=':8,'Ndn=':8,'D=':8);

Nup:=0; Ndn:=0;

for p:=1 to ps do

begin

nc[p,1]:=0; nc[p,2]:=0; d[p]:=0;

for kR:=1 to 2*ks do

begin

if kR<=ks/2 then nc[p,1]:=nc[p,1]+nk1[kR,p];

if (kR>ks/2) and (kR<=1.5*ks) then nc[p,2]:=nc[p,2]+nk1[kR,p];

if kR>1.5*ks then nc[p,1]:=nc[p,1]+nk1[kR,p];

if kR<=ks then d[p]:=d[p]+nk1[kR,p];

if kR>ks then d[p]:=d[p]-nk1[kR,p];

end;

d[p]:=d[p]/n/q1/ks;

nc[p,1]:=nc[p,1]/n/q1/ks; nc[p,2]:=nc[p,2]/n/q1/ks;

write(out2,nc[p,1]:8:4,nc[p,2]:8:4,d[p]:8:4); writeln(out2);

Nup:=Nup+nc[p,1]; Ndn:=Ndn+nc[p,2];

end;

{ vvod(Nsu);}

Nup:=Nup/ps; Ndn:=Ndn/ps;

writeln(out2,'Conc':8,'priv':8,'volms':8);

for p:=1 to ps do

begin

write(out2,(nk1[1,p]+nk1[2*ks,p])/2/n/q1:8:3);

for kR:=2 to 2*ks do write(out2,(nk1[kR,p]+nk1[kR-1,p])/2/n/q1:8:3);

writeln(out2);

end;

summa;

writeln(out2);

close(out2);

Assign(out3,'c:\80int.txt'); Append(out3);

write(out3,ww0:9:4,vvsr/n:9:4,vsr/n:9:4,Tsr/n:9:4,Tss:9:4,Fdir:9:4,

conc:9:4,NNN:9:4,Ndn:9:4,Nup:9:4,Ndn-Nup:9:4);

write(out3,nc[2,2]-nc[2,1]:9:4); {changeble addon}

for p:=1 to 2 do

for kL:=1 to ps do write(out3,nc[kL,p]:9:4);

writeln(out3);

close(out3);

until m1=0;

writeln('Ready'); close(inp1);

end.

ПРИЛОЖЕНИЕ 2

Дата добавления: 2014-11-29