unit Unit_cilindr;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TPoint3D = Record //Точка в трехмерном пространстве
   X, Y, Z : Integer;
  End;

  TfmForm1 = class(TForm)
    Button1: TButton;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  fmForm1: TfmForm1;

  Arr : Array[1..50, 1..50, 1..50] Of Integer; //Массив из 50х50х50 элементов
  Napr1 : Array[1..50, 1..50, 1..50] Of TPoint3D;
        //50х50х50 (для каждой точки) направляющих векторов длиной = 1 - для
        // поиска координат следующей точки
  E : Array[1..50, 1..50, 1..50] Of Double; //Напряженность для каждой точки

implementation

{$R *.dfm}

procedure TfmForm1.Button1Click(Sender: TObject);
Var
 Cycle, Cycle2, Cycle3,k : Integer; //Для организации циклов
 LengthOfNapr : Double; //Длина направляющего вектора Napr
 H,dh,r,S,ds : Double; //Длина вектора Napr1, так как она не всегда равна 1
 Q1,Q,C : Double; //Общий заряд при dS = h*h = 1*1
 Q2,C1,C2 : Double; //Общий заряд при dS = H*H
begin
 For Cycle := 1 To 50 Do
  For Cycle2 := 1 To 50 Do
   For Cycle3 := 1 To 50 Do
    Arr[Cycle, Cycle2, Cycle3] := 0;
 k:=0;
 S:=0;
 ds:=0;
 r:= 0.0013;
 For Cycle := 1 To 50 Do
  For Cycle2 := 1 To 50 Do
   For Cycle3 := 1 To 50 Do
   Begin
    If ((Cycle-25) * (Cycle-25) + (Cycle2-25) * (Cycle2-25) + (Cycle3-25) * (Cycle3-25) <= 169) And
       ((Cycle-25) * (Cycle-25) + (Cycle2-25) * (Cycle2-25) + (Cycle3-25) * (Cycle3-25) > 144)Then
             //Проверка на принадлежность сфере - все точки в пределах (12, 13]
     begin
      k:=k+1;
      Arr[Cycle, Cycle2, Cycle3] := 1;
     end;        //Устанавливаем в 1 - как признак того, что точка принадлежит сфере
   End;
 S:=4*pi*r*r;
 ds:=S/k;
 dh:=0.0001;
 Q1 := 0;
 Q2 := 0;
 C1 := 0;
 C2 := 0;
 For Cycle := 1 To 50 Do
  For Cycle2 := 1 To 50 Do
   For Cycle3 := 1 To 50 Do
    If Arr[Cycle, Cycle2, Cycle3] = 1 Then
    Begin
     LengthOfNapr := Sqrt(Sqr(Cycle - 25) + Sqr(Cycle2 - 25) + Sqr(Cycle3 - 25));
                     //Корень из суммы квадратов: Sqr - квадрат, Sqrt - корень
     //Координаты точек направляющего вектора:
     // начало - центр сферы (25, 25, 25)
     // конец - точка на сфере (Cycle, Cycle2, Cycle3)

     Napr1[Cycle, Cycle2, Cycle3].X := Round((Cycle - 25) / LengthOfNapr);
     Napr1[Cycle, Cycle2, Cycle3].Y := Round((Cycle2 - 25) / LengthOfNapr);
     Napr1[Cycle, Cycle2, Cycle3].Z := Round((Cycle3 - 25) / LengthOfNapr);

     H := Sqrt(Sqr(Napr1[Cycle, Cycle2, Cycle3].X) +
               Sqr(Napr1[Cycle, Cycle2, Cycle3].Y) +
               Sqr(Napr1[Cycle, Cycle2, Cycle3].Z))*dh;

     If H <> 0 Then
     Begin
      E[Cycle, Cycle2, Cycle3] := 2 / H;
      Q1 := Q1 + E[Cycle, Cycle2, Cycle3] * ds;
      Q2 := Q2 + E[Cycle, Cycle2, Cycle3]* dh * dh;
     End;
    End;
     { Q2:=(Q2*24.6176e-6)/dh*dh;}
      Q1:=Q1*8.854e-12;
      Q2:=Q2*8.854e-12;
      C1:=Q1/20;
      C2:=Q2/20;
      C:=4*pi*r*8.854e-12*1.00058;
 ShowMessage('       Q1 = ' + FormatFloat('0.000E+00     ', Q1) + #10#13 +
             '       Q2 = ' + FormatFloat('0.000E+00     ', Q2)+#10#13 +'По теореме Гаусса (ds=h*h):   C1 = ' +FormatFloat('0.000E+00', C1) + #10#13 +
             'По теореме Гаусса (ds=S/n):   C2 = ' + FormatFloat('0.000E+00     ', C2)+ #10#13 +
             'По формуле (С=4*pi*eps*R):     C   = ' + FormatFloat('0.000E+00      ', C));
 end;

end.