Modmd82krug

 

program dat(input,output);

 label ll; const ks=12; ps=3; ProgID=82; r21=0.600;r22=0.240;

 var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1: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: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>r21-r22)and(r0<r21+r22) then

  begin

    sum:=(r21+r22+r0)/2;

    y0:=4/pi*2*arctan(sqrt((sum-r22)*(sum-r21)*(sum-r0)/sum)/(sum-r22));

    y:=random;

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

  end;

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*r22*r22)*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);

            {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)

              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;

     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;

     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.

ПРИЛОЖЕНИЕ 3