program field; {$n+,e-,m 65520,260000,260000} uses crt,dos; type Uz=array[1..7,0..4] of extended; Werte=record hf,s,y,u0,i,j,k,o,w,z,i1,i2,i3,i4: integer; dg0,z0,dr,r,r0,dt,t,su,dg,K0,c1,c2,c3,c4: extended; c01,c02,c03,c04,c10,c20,c30,c40,c1s,c2s,c3s,c4s: extended; U: array[0..2] of Uz; Ug,Uh: array[1..7,1..730] of single; z1: single; ul: array[1..4] of extended; V: array[1..7] of extended; a,q: array[1..4] of extended; g0,h0,m0,g,h,m,p,f,ri: array[1..4,1..4] of extended; n0,l,n: array[1..4,1..4,1..4] of extended; Cs,Tg: array[0..730] of extended; x1: char; a0: double; xi,code: integer; bn,bs: Boolean end; var wr: Werte; ca: char; x2:string[3]; fi1: text; fi6:file of extended; fi2: file of double; fi3:file of single; fi5: file of Werte; bo: Boolean; Us: array[0..730] of ^Uz; fi4: file of Uz; ExitSave: Pointer; {$F+} procedure MyExit; {$F-} begin close(fi1); if wr.bn then begin close(fi2); close(fi3); if wr.bs then close (fi6) end; ClrScr; ExitProc:= ExitSave; end; procedure kSchleife; begin with wr do begin for i:=1 to 4 do if (i=1) or (i=4) then begin g0[i,i]:=U[1][i,2]; l[i,i,1]:=(U[1][i,1]-U[1][i,3])/2/dr; l[i,i,2]:=(U[2][i,2]-U[0][i,2])/2/dt; m0[i,i]:=(U[2][i,2]-2*U[1][i,2]+U[0][i,2])/sqr(dt); g[i,i]:=g0[i,i]+1;m[i,i]:=m0[i,i] end; g0[2,2]:=r*r*U[1][2,2]; l[2,2,1]:=2*r*U[1][2,2]+r*r*(U[1][2,1]-U[1][2,3]) /2/dr; l[2,2,2]:=r*r*(U[2][2,2]-U[0][2,2])/2/dt; m0[2,2]:=r*r*(U[2][2,2]-2*U[1][2,2]+U[0][2,2])/sqr(dt); g0[3,3]:=sqr(r*Cs[k])*U[1][3,2]; l[3,3,1]:=sqr(r*Cs[k])*(2*U[1][3,2]/r +(U[1][3,1]-U[1][3,3])/2/dr); l[3,3,2]:=sqr(r*Cs[k])*((U[2][3,2]- U[0][3,2])/2/dt-2*U[1][3,2]*Tg[k]); g0[3,4]:=U[1][5,2]*sqr(Cs[k])/r; g0[4,3]:=g0[3,4]; a[3]:=U[1][6,2]*sqr(Cs[k])/r; a[4]:=U[1][7,2]/r; l[3,4,1]:=sqr(Cs[k])*(-U[1][5,2]/r+(U[1][5,1]-U[1][5,3])/2/dr)/r; l[3,4,2]:=sqr(Cs[k])*((U[2][5,2]-U[0][5,2])/2/dt-2*U[1][5,2]* Tg[k])/r; l[4,3,1]:=l[3,4,1]; l[4,3,2]:=l[3,4,2]; p[3,1]:=sqr(Cs[k])*(-U[1][6,2]/r+(U[1][6,1]-U[1][6,3])/2/dr)/r; p[3,2]:=sqr(Cs[k])*((U[2][6,2]-U[0][6,2])/2/dt-2*U[1][6,2]*Tg[k]) /r; p[4,1]:=-U[1][7,2]/r/r+(U[1][7,1]-U[1][7,3])/2/dr/r; p[4,2]:=(U[2][7,2]-U[0][7,2])/2/dt/r; m0[3,4]:=sqr(Cs[k])*((U[2][5,2]-2*U[1][5,2]+U[0][5,2])/dt/dt -2*(U[2][5,2]-U[0][5,2])/dt*Tg[k]+2*U[1][5,2]*(sqr(Tg[k])-1))/r; m0[4,3]:=m0[3,4]; q[3]:=sqr(Cs[k])*((U[2][6,2]-2*U[1][6,2]+U[0][6,2])/sqr(dt) -2*(U[2][6,2]-U[0][6,2])/dt*Tg[k]+2*U[1][6,2]*(sqr(Tg[k])-1))/r; q[4]:=(U[2][7,2]-2*U[1][7,2]+U[0][7,2])/sqr(dt)/r; g[2,2]:=g0[2,2]+r*r; m[2,2]:=m0[2,2]; g[3,3]:=g0[3,3]+sqr(r*Cs[k]); g[3,4]:=g0[3,4];g[4,3]:=g0[3,4]; m[3,4]:=m0[3,4];m[4,3]:=m0[3,4] {Potentiale und Ableitungen}; if (g[1,1]<-100) or (g[2,2]<-100*r*r) or (g[3,3]<-100*sqr(r*Cs[k])) or (g[4,4]<-100) then hf:=y; if (abs(g[3,4])>100*r*Cs[k]) or (abs(a[3])>100*r*Cs[k]) or (abs(a[4])>100) then hf:=y; dg0:=g0[3,3]+sqr(r*Cs[k])*g0[4,4]+g0[3,3]*g0[4,4]+sqr(g0[3,4]); dg:=dg0+sqr(r*Cs[k]); h0[1,1]:=(r*r*dg0+g0[2,2]*dg)/sqr(r*r*Cs[k]); h[1,1]:=h0[1,1]+1; h0[2,2]:=(dg0+g0[1,1]*dg)/sqr(r*r*Cs[k]); h[2,2]:=h0[2,2]+1/r/r; h0[3,3]:=(g[1,1]*g[2,2]*g0[4,4]+g0[1,1]*g[2,2]+g0[2,2]) /sqr(r*r*Cs[k]);h[3,3]:=h0[3,3]+1/sqr(r*Cs[k]); h0[4,4]:=(g0[1,1]*g[2,2]*g[3,3]+g0[2,2]*g[3,3]+r*r*g0[3,3]) /sqr(r*r*Cs[k]);h[4,4]:=h0[4,4]+1; h0[3,4]:=-g[1,1]*g[2,2]*g[3,4]/sqr(r*r*Cs[k]); h0[4,3]:=h0[3,4]; h[3,4]:=h0[3,4]; h[4,3]:=h[3,4]; for i:=1 to 4 do for j:=1 to 4 do f[i,j]:=p[i,j]-p[j,i]; su:=0; for i:=1 to 4 do for j:=1 to 4 do for o:=1 to 4 do for w:=1 to 4 do begin if (i=4) or (j=4) or (o=4) or (w=4) then s:=-1 else s:=1; su:=su+s*h[i,j]*h[o,w]*f[i,o]*f[j,w] end; for i:=1 to 4 do for j:=1 to 4 do begin ri[i,j]:=g[i,j]*(su/4-3*K0); for o:=1 to 4 do for w:=1 to 4 do begin if (i=4) and (j=4) or (o=4) or (w=4) then s:=-1 else s:=1; ri[i,j]:=ri[i,j]-s*h[o,w]*f[i,o]*f[j,w] end end; for i:=1 to 4 do for j:=1 to 4 do for o:=1 to 4 do begin n0[i,j,o]:=0; for w:=1 to 4 do begin if (o=3) and (w=4) and ((i=3) or (j=3)) or (o=4) and (w=3) and ((i=4) or (j=4)) then s:=-1 else s:=1; n0[i,j,o]:=n0[i,j,o]+s* h[o,w]*(l[j,w,i]+l[i,w,j]-l[i,j,w])/2 end; n[i,j,o]:=n0[i,j,o] end; n[2,2,1]:=n0[2,2,1]-r*h[1,1]; n[1,2,2]:=n0[1,2,2]+r*h[2,2]; n[2,1,2]:=n[1,2,2]; n[3,3,1]:=n0[3,3,1]-r*sqr(Cs[k])*h[1,1]; n[1,3,3]:=n0[1,3,3]+r*sqr(Cs[k])*h[3,3]; n[3,1,3]:=n[1,3,3]; n[3,3,2]:=n0[3,3,2]+Tg[k]*sqr(r*Cs[k])*h[2,2]; n[2,3,3]:=n0[2,3,3]-Tg[k]*sqr(r*Cs[k])*h[3,3]; n[3,2,3]:=n[2,3,3]; V[1]:=2*n[1,1,1]*(n[1,1,1]-1/r)+n[1,1,2]*(Tg[k]+ 2*n[2,2,2])+m[1,1]*h[2,2]/2-ri[1,1]; V[1]:=V[1]+r*h0[2,2]*n[1,2,2]+n0[1,2,2]/r+h0[2,2]; V[1]:=V[1]+r*sqr(Cs[k])*h0[3,3]*n[1,3,3] +n0[1,3,3]/r+sqr(Cs[k])*h0[3,3]; for i:= 1 to 4 do for j:= 1 to 4 do begin if (i=3) and (j=4) or (i=4) and (j=3) then s:=-1 else s:=1; V[1]:=V[1]+s*n0[1,i,j]*n[1,j,i] end; V[1]:=2*g[1,1]*V[1]; V[2]:=2*(-n[2,2,1]*n[1,1,1]+n0[2,2,1]/r -h0[1,1])+ri[2,2]; V[2]:=V[2]-n[2,2,2]*(2*n[2,2,2]+Tg[k])+h[2,2]*m[2,2]/2 -r*h0[2,2]*n[2,2,1]+r*h0[1,1]*n[2,1,2]; V[2]:=V[2]+Tg[k]*sqr(r*Cs[k])*h0[3,3]*n[2,3,3] -n0[2,2,1]/r+h0[1,1]; V[2]:=V[2]+r*(n0[2,1,2]+r*h0[2,2])+Tg[k]*(n0[2,3,3] -Tg[k]*sqr(r*Cs[k])*h0[3,3]); for i:=1 to 4 do for j:=1 to 4 do begin if (i=3) and (j=4) or (i=4) and (j=3) then s:=-1 else s:=1; V[2]:=V[2]-s*n0[2,i,j]*n[2,j,i] end; V[2]:=2*g0[1,1]+2*g[1,1]*V[2]; V[4]:=2*g[1,1]*(2*n[4,4,1]*(1/r-n[1,1,1])-n[4,4,2]*(Tg[k]+ 2*n[2,2,2])-m[4,4]*h[2,2]/2+ri[4,4]); for i:=1 to 4 do for j:= 1 to 4 do begin if (i=3) or (j=3) then s:=-1 else s:=1; V[4]:=V[4]-s*2*g[1,1]*n[4,i,j]*n[4,j,i] end; V[5]:=2*g[1,1]*(2*n[3,4,1]*(1/r-n[1,1,1])-n[3,4,2]*(Tg[k]+2*n[2,2,2] )-m[3,4]*h[2,2]/2+ri[3,4]); for i:=1 to 4 do for j:=1 to 4 do V[5]:=V[5]-2*g[1,1]*n[3,i,j]*n[4,j,i]; for w:=6 to 7 do begin V[w]:=g[1,1]*(2*n[1,1,1]*h[1,1]*f[w-3,1]+2*n[2,2,2]*h[2,2]* f[w-3,2]-q[w-3]*h[2,2]); for i:=1 to 4 do for j:=1 to 4 do for o:=1 to 4 do begin if (w=6) and ((i=4) or (j=4) or (o=4)) then s:=-1 else s:=1; V[w]:=V[w]+g[1,1]*(s*n[w-3,i,j]*f[j,o] -n[i,j,j]*f[w-3,o])*h[o,i] end end; for i:=1 to 4 do if (i=1) or (i=4) then U[1][i,4]:=2*U[1][i,2]-U[1][i,0]+V[i]*4*dr*dr; U[1][2,4]:=2*U[1][2,2]-U[1][2,0]+sqr(2*dr/r)*(V[2]- 2*U[1][2,2]-2*r*(U[1][2,1]-U[1][2,3])/dr); U[1][7,4]:=2*U[1][7,2]-U[1][7,0]+sqr(2*dr)*V[7]*r- 8*U[1][7,2]*sqr(dr/r)+4*(U[1][7,1]-U[1][7,3])/dr/r; for i:=5 to 6 do U[1][i,4]:=2*U[1][i,2]-U[1][i,0]+(U[1][i,1] -U[1][i,3])*4*dr/r-U[1][i,2]*8*sqr(dr/r)+V[i]*sqr(2*dr/Cs[k])*r; U[1][3,4]:=-((1+U[1][1,4])*(1+U[1][2,4])*(U[1][4,4]+sqr(U[1][5,4] *Cs[k]/r/r))+U[1][2,4]*(1+U[1][1,4])+U[1][1,4]); U[1][3,4]:=U[1][3,4]/((1+U[1][4,4])*(1+U[1][1,4])*(1+U[1][2,4])) end; end; begin { main program } ExitSave := ExitProc; ExitProc := @MyExit; with wr do begin for k:=0 to 730 do GetMem(Us[k],350); write(#13#10'file Nr. (Anfang): '); ReadLn(x2); write('Suchen (U), Abgebrochenes weiter (W), Speicher (S): '); ca:=ReadKey; WriteLn(ca); if upcase(ca)='W' then bo:=false else bo:=true; if upcase(ca)='U' then bn:=false else bn:=true; if upcase(ca)='S' then bs:=true else bs:=false; if not bo then begin assign(fi5,'c:\data\fsp'+x2+'.dat'); reset(fi5); read(fi5,wr); close(fi5); erase(fi5); assign(fi4,'c:\data\fz'+x2+'.dat'); reset(fi4); for i:=0 to z do read(fi4,Us[i]^); close(fi4); erase(fi4); assign(fi1,'c:\data\fet'+x2+'.dat'); append(fi1); if bn then begin if bs then begin assign(fi6,'c:\data\fs'+x2+'.dat'); reset(fi6); seek(fi6,FileSize(fi6)) end; assign(fi2,'c:\data\fe'+x2+'.dat'); reset(fi2); seek(fi2,FileSize(fi2)); assign(fi3,'c:\data\fgr'+x2+'.dat'); reset(fi3); seek(fi3,FileSize(fi3)) end end else begin write('Feinheit (z<=730) z='); ReadLn(z); z0:=z/1; write('Anfangsradius r0='); ReadLn(r0); write('Schrittweite dr='); ReadLn(dr); write('konst.Kruemmung K0='); ReadLn(K0); if bn then begin write('Automatik (a): '); x1:=ReadKey; WriteLn(x1) end else x1:='A'; if upcase(x1)='A' then begin write('"Kleinstmasse" m min = '); ReadLn(c01); write('"Groesstmasse" m max = '); ReadLn(c10); write('step '); ReadLn(c1s); write('"Kleinstdrehimpuls" s min = '); ReadLn(c02); write('"Groesstdrehimpuls" s max = '); ReadLn(c20); write('step '); ReadLn(c2s); write('"Kleinstladung" Q min = '); ReadLn(c03); write('"Groesstladung" Q max = '); ReadLn(c30); write('step '); ReadLn(c3s); write('"kleinst.magn.Moment" M min = '); ReadLn(c04); write('"groesst.magn.Moment" M max = '); ReadLn(c40); write('step '); ReadLn(c4s); c1:=c01;c2:=c02;c3:=c03;c4:=c04 end else begin write('"Masse" m='); ReadLn(c1); write('"Drehimpuls" s='); ReadLn(c2); write('"Ladung" Q='); ReadLn(c3); write('"magn.Moment" M='); ReadLn(c4) end; write('Abbruch mit ^U.'); delay(2000); k:=0; dt:=pi/(2*z-3); repeat t:=(k-1)*dt; Cs[k]:=cos(t); Tg[k]:=sin(t)/Cs[k]; inc(k); until k=z+1; for i:=1 to 4 do begin a[i]:=0; q[i]:=0; for j:=1 to 4 do begin g[i,j]:=0; h[i,j]:=0; p[i,j]:=0; m[i,j]:=0; g0[i,j]:=0; h0[i,j]:=0; m0[i,j]:=0; for o:=1 to 4 do l[i,j,o]:=0 end end; z1:=(z-1)/1 end; repeat clrscr; writeln('z=',z:2,',r0=',r0:3:1,',dr=',dr:5:3,',K0=',K0:8,','); writeln('m=',c1:11,',s=',c2:11,',Q=',c3:11,',M=',c4:11); window(1,3,80,25); if bo then begin hf:=0; assign(fi1,'c:\data\fet'+x2+'.dat'); if bn then begin if bs then begin assign(fi6,'c:\data\fs'+x2+'.dat'); rewrite(fi6) end; rewrite(fi1); assign(fi2,'c:\data\fe'+x2+'.dat'); rewrite(fi2); assign(fi3,'c:\data\fgr'+x2+'.dat'); rewrite(fi3); write(fi3,z1) end else if (c1=c01) and (c2=c02) and (c3=c03) and (c4=c04) then rewrite(fi1) else append(fi1); writeln(fi1,'z=',z:4,',r0=',r0:5:3,',dr=',dr:6:4,',K0=',K0:8,','); write(fi1,'m=',c1:12,',s=',c2:12,',Q=',c3:12,',M=',c4:12,','); k:=0; repeat y:=0; repeat r:=r0-y*dr; Us[k]^[1,y]:=K0*r*r+c1/r-sqr(c3/r)/2+sqr(c4/r/r)*(1+sqr(Cs[k]))/10; Us[k]^[2,y]:=sqr(c4/r/r)*(sqr(Cs[k])/3-0.3); Us[k]^[3,y]:=sqr(c4/r/r)*(sqr(Cs[k])/15-0.3); Us[k]^[4,y]:=-K0*r*r-c1/r+sqr(c3/r)/2+sqr(c4/r/r)*(1-sqr(Cs[k]))/2; Us[k]^[5,y]:=c2-c3*c4/r/2; Us[k]^[6,y]:=c4; Us[k]^[7,y]:=c3; inc(y); until y=4; inc(k) until k=z+1 {Anfangswerte gesetzt}; for y:=0 to 3 do begin r:=r0-y*dr; for i:=1 to 4 do begin write(r:6:3,' ç',i,i); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; a0:=Us[k]^[i,y]; write(a0:8); if bn then write(fi2,a0) end end; write(r:6:3,' ç34'); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; a0:=Us[k]^[5,y]*Cs[k]/r/r; write(a0:8); if bn then write(fi2,a0) end; write(r:6:3,' 3'); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; a0:=Us[k]^[6,y]*Cs[k]/r/r; write(a0:8); if bn then write(fi2,a0) end; write(r:6:3,' A4'); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; a0:=Us[k]^[7,y]/r; write(a0:8); if bn then write(fi2,a0) end end; for y:=1 to 7 do for k:=1 to 730 do begin Uh[y,k]:=0; Ug[y,k]:=0 end; y:=2 end; repeat if bo then begin r:=r0-y*dr; k:=1 end; repeat bo:=true; if KeyPressed then if ReadKey=#21 then begin assign(fi5,'c:\data\fsp'+x2+'.dat'); rewrite(fi5); write(fi5,wr); close(fi5); assign(fi4,'c:\data\fz'+x2+'.dat'); rewrite(fi4); for i:=0 to z do write(fi4,Us[i]^); close(fi4); halt end; for o:=0 to 2 do U[o]:=Us[k+o-1]^; kSchleife; Us[k]^:=U[1]; inc(k) until k=z; for i:=1 to 7 do begin Us[0]^[i,4]:=Us[2]^[i,4]; Us[z]^[i,4]:=Us[z-1]^[i,4] end {"Randwerte" geladen}; for k:=0 to z do begin U[1]:=Us[k]^; for i:=1 to 7 do begin if bs and bn then write(fi6,U[1][i,4]); for j:=0 to 3 do U[1][i,j]:=U[1][i,j+1] end; Us[k]^:=U[1] end {nachgerueckt}; r:=r-2*dr; for k:=1 to z-1 do begin U[1]:=Us[k]^; ul[1]:=K0*r*r+c1/r-sqr(c3/r)/2+sqr(c4/r/r)*(1+sqr(Cs[k]))/10; ul[2]:=sqr(c4/r/r)*(sqr(Cs[k])/3-0.3); ul[3]:=sqr(c4/r/r)*(sqr(Cs[k])/15-0.3); ul[4]:=-K0*r*r-c1/r+sqr(c3/r)/2+sqr(c4/r/r)*(1-sqr(Cs[k]))/2; for i:=1 to 4 do begin if abs(U[1][i,4])>1e12 then hf:=y; if (Uh[i,k]=0) and (U[1][i,4]<-0.9) then Uh[i,k]:=r; if (Ug[i,k]=0) and (abs(U[1][i,4])>2*abs(ul[i])) and (ul[i]<>0) then Ug[i,k]:=r end; for i:=5 to 6 do begin if abs(U[1][i,4])*Cs[k]>1e12*r*r then hf:=y; if (Uh[i,k]=0) and (abs(U[1][i,4])*Cs[k]>0.9*r*r) then Uh[i,k]:=r end; if (Ug[5,k]=0) and (abs(U[1][5,4])>2*abs(c2-c3*c4/r/2)) and (c2-c3*c4/r/2<>0) then Ug[5,k]:=r; if (Ug[6,k]=0) and (abs(U[1][6,4])>2*abs(c4)) and (c4<>0) then Ug[6,k]:=r; if (Ug[7,k]=0) and (abs(U[1][7,4])>2*abs(c3)) and (c3<>0) then Ug[7,k]:=r; if (Uh[7,k]=0) and (abs(U[1][7,4])>0.9*r) then Uh[7,k]:=r; if abs(U[1][7,4])>1e12*r then hf:=y end; for i:=1 to 4 do begin write(r:6:3,' ç',i,i); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; a0:=Us[k]^[i,4]; write(a0:8); if bn then write(fi2,a0) end end; write(r:6:3,' ç34'); for u0:=0 to 6 do begin k:=round(u0*(z-1.5)/6)+1; if abs(r)0 then y:=round(r0/dr) until y=round(r0/dr); if hf=0 then writeln(fi1,y+1:4) else writeln(fi1,hf+2:4); close(fi1); if bn then begin for y:=1 to 7 do begin for k:=1 to z-1 do write(fi3,Ug[y,k]); for k:=1 to z-1 do write(fi3,Uh[y,k]) end; if bs then begin write(fi6,z0,z0,z0,z0,z0,z0,z0); close(fi6) end; close(fi2); close(fi3) end; if upcase(x1)='A' then begin c1:=c1+c1s; if c1>c10 then begin c1:=c01; c2:=c2+c2s end; if c2>c20 then begin c2:=c02;c3:=c3+c3s end; if c3>c30 then begin c3:=c03;c4:=c4+c4s end; if bn then begin val(x2,xi,code); inc(xi); str(xi,x2) end else code:=0 end; window(1,1,80,25) until (c4>c40) or (upcase(x1)<>'A') or (code<>0); for k:=0 to 730 do FreeMem(Us[k],350); reset(fi1); if bn then begin reset(fi2); reset(fi3); if bs then reset(fi6) end end; end.