Uses Crt, Math, Printer; Const N1=1200; N=2; // g: Double = 9.81; //przyspieszenie ziemskie [m/sek2] po: Double = 1.0; //cisnienie odniesienia [kGm/cm2] eps: Double = 1.0E-15; hmin: Double =1E-15; Type Wek= Array[0..N1,0..3] of Double; Wektor= Array[0..N] of Double; Var x, xe, h:Double; fi: Boolean; Y: Wektor; A: Wek; i, j: Longint; k, l, m: Integer; // Col, Lin: Integer; { Procedura zapisuje w pliku "Nazwa.*" wyniki obliczeń z programu SCRATONUR:} Procedure DyskWrit (N: Integer; A: Wek); Var Nazwa: String[36]; Plik: File of Double; Begin // Readln; Write('Podaj nazwe zbioru na dysku = '); Readln(Nazwa); Nazwa:='C:\LAZARUS\PAPIEROS\Dane\'+Nazwa; Assign(Plik,Nazwa); Reset(Plik); A[0,0]:= N+0.0; Write(Plik,A[0,0],A[0,1],A[0,2]); For k:=1 To N Do For l:=0 To 2 Do Write(Plik,A[k,l]); Close(Plik) End; //Dyskwrit {Program odczytuje z pliku "Nazwa.*" wyniki obliczeń z programu SCRATONUR: N= liczba pomiarów, x= z, Y[1] = p, Y[2] = t} Procedure DyskRead; Var k, l: Integer; N: Integer; Nazwa: String[30]; Plik: File of Double; Begin Write('Podaj nazwe zbioru na dysku = '); Readln(Nazwa); Nazwa:='C:\LAZARUS\PAPIEROS\Dane\'+Nazwa; Assign(Plik,Nazwa); Rewrite(Plik); Read(Plik, A[0,0]); N:= Round(A[0,0]); For k:=1 To N Do For l:=0 To 2 Do Read(Plik,A[k,l]); Close(Plik); ClrScr; // Writeln(#10,'Liczba pomiarów N= ',N,#10,#13); For i:=0 To N Do Writeln('z= ',A[i,0],' p=',A[i,1],' t=',A[i,1]); End; Procedure Fun(N:Byte; z:Double; Y:Wektor; Var K:Wektor); Const R: Double = 1.0; //PROMIEŃ KOMORY SPALANIA [cm] L: Double = 10.0; //DLUGOSC KOMORY SPALANIA [cm] Rkr: Double = 0.25; //PROMIEŃ KRYTYCZNY DYSZY [cm] Ro: Double = 1.3; //GĘSTOSC PALIWA [g/cm3] fo: Double = 0.23451E3; //SILA PALIWA [kGcm/g] Uo: Double = 1.72; //LINIOWA SZYBKOSC SPALANIA [cm/sek] np: Double = 0.2;//0.172; //WYKLADNIK RÓWNANIA SPALANIA 0.164+/-0.017 ka: Double = 1.2; //WYKLADNIK ADIABATY Var S, Skr, Vo ,Vk: Double; a, b, e, eo: Double; xk, Delta, C, mp: Double; pgr, M1, M2, M3, M4: Double; p1, p2: Double; Function Zs(z1, z2, z: Double): Double; // Funkcja Zs dla z2 z2 Then Zs := 0.0 End; //Zs Function Zp(z1, z2, z: Double): Double; // Funkcja Zp dla z2 z2 Then Zp := z2 Else If z > 1 Then Zp := 0 End; //Zp Function f(p: Double): Double; Begin f:= 1.0- M2*p**(1.0-np)- M3*p End; Begin // If fi = False Then Begin eo:= L; a:= eo/R; b:= eo/Rkr; e:= eo/L; Vo:= Pi*eo*eo*eo/(a*a*e); Vk:= Pi*eo*eo*eo*(1.0/(a*a*e)+(a*a+b*b+a*b)*(b-a)/(3*a*a*a*b*b*b)); S:= Pi*eo*eo/(a*a); Skr:= Pi*eo*eo/(b*b); C:= 0.95(**Sqrt(g)*)*(2.0/(ka+1.0))**(1.0/(1.0-ka)) * Sqrt(2.0*ka/(ka+1.0))/Sqrt(fo); // Psi:= z; mp:= Ro*Vo; Delta:= mp/Vk; xk:= Ro/Delta; M1:= fo*Ro/po; M2:= C*Skr/(Ro*(S*Zs(0.0,1.0,z))*Uo); M3:= 1.0/(Ro*fo); M4:= mp/(Ro*(S*Zs(0.0,1.0,z))*Uo); If fi=False Then Begin p1:= 0.0; p2:= (1.0/M2)**(1.0/(1.0-np)); Repeat pgr:= (p1+p2)/2.0; If f(pgr)*f(p2)<0.0 Then p1:= pgr Else p2:= pgr; Until (Abs(p1-p2) < 1E-10); End; K[1]:= M1*(1.0- M2*(Y[1]/po)**(1.0-np)- M3*(Y[1]/po))/(Zp(0.0,1.0,z)+xk); K[2]:= M4/((Y[1]/po)**np); // GotoXY(50,1); Writeln(' pgr= ',pgr:4:6,' Ciag= ',1.5*Skr*(pgr/po):4:6); End; Procedure ScratonUR(N:Byte; Var x:Double; xe:Double; Var Y:Wektor); Label A,B,C; Var i: Byte; E,K1,K2,K3,K4,K5,W: Wektor; E1,K2M: Double; Begin W:=Y; If fi Then Begin A: E1:=Abs((xe-x)/h); If E1<=1.5 Then fi:=E1>1.0 Else Goto C End; If fi Then h:=0.5*(xe-x) Else h:=xe-x; C: Fun(N,x,W,K1); For i:=1 To N Do W[i]:=Y[i]+2.0*K1[i]/9.0; Fun(N,x+2.0*h/9.0,W,K2); For i:=1 To N Do W[i]:=Y[i]+K1[i]/12.0+0.25*K2[i]; Fun(N,x+h/3.0,W,K3); For i:=1 To N Do W[i]:=Y[i]+0.5390625*K1[i]-1.8984375*K2[i]+2.109375*K3[i]; Fun(N,x+0.75*h,W,K4); For i:=1 To N Do W[i]:=Y[i]-0.3105*K1[i]+1.8225*K2[i]-1.1016*K3[i]+0.4896*K4[i]; Fun(N,x+0.9*h,W,K5); K2M:=0.0; For i:=1 To N Do If K4[i]<>K1[i] Then Begin E[i]:=(0.34*K1[i]-0.972*K2[i]+1.632*K4[i]-K5[i])*(1.29357298475E-1*K1[i]-5.51470588235E-1*K2[i]+0.46568627451*K3[i]-4.35729847494E-2*K4[i])/(K4[i]-K1[i]); E1:=Abs(E[i]); If E1>K2M Then K2M:=E1 End Else E[i]:=0.0; // Writeln(' E[1]=',E[1] ,' E[2]=',E[2],' K2M=',K2M); // If K2M>eps Then Goto B; Inc(j); For i:=1 To N Do Y[i]:=Y[i]-E[i]+1.04938271605E-1*K1[i]+4.76470588235E-1*K3[i]+2.37037037037037E-1*K4[i]+1.81554103123E-1*K5[i]; If fi Then Begin x:=x+h; B: If K2M<>0.0 Then h:= h*Exp(0.2*Ln(eps/K2M)) Else h:= hmin/2.0; If Abs(h)>hmin Then Goto A; // (*ClrScr; GotoXY(10,10);*) Writeln('Nie mozna uzyskac dokladnosci eps= ',eps); // Repeat Until KeyPressed; h:=h/10.0; Goto C End; x:=xe; fi:=True End; {ScratonUR} Begin // Warunki początkowe x,y fi:= False; Y[1]:= 10.0; Y[2]:= 0.0; x:= 0.0; xe:= x+0.001; i:= 0; j:= 0; m:= 1; A[0,0]:= 1001.0; A[0,1]:= 0.0; A[0,2]:= 0.0; Writeln(' x',' t',' p',' j',#10); // Col:= WhereX; Lin:= WhereY; While x<= 1.001 Do Begin A[m,0]:= x; A[m,1]:= Y[2]; A[m,2]:= Y[2]; DyskWrit(1002,A); If i <= 25 then Begin Writeln(' ',x:4:3,' ',Y[2]:4:6 ,' ',Y[1]:4:6,' ',j); Inc(i); Inc(m); If i > 25 then Begin (*Delay(10000)*); i:= 0; ClrScr; Writeln(' x',' t',' p',' j',#10); End; End; ScratonUR(N,x,xe,Y); xe:=x+0.001; End; DyskRead; Repeat Until KeyPressed; End.