勒让德-高斯求积法求磁感应强度.pas

勒让德-高斯求积法求磁感应强度

unit Unit24;

interface
uses
  Math;

  function FLRGS(A,B:extended;F1,F2:TCmplxVector;Frequency:array of extended):extended;
  //勒让德-高斯求积法

implementation

function FLRGS(A,B:extended;F1,F2:TCmplxVector;Frequency:array of extended):extended;
//勒让德-高斯求积法
const
  IntegPoint:array[1..7]of extended=(-0.9491079123,-0.7415311856,-0.4058451514,0.0,0.4058451514,0.7415311856,0.9491079123);
  coef:array[1..7]of extended=(0.1294849662,0.2797053915,0.3818300505,0.4179591837,0.3818300505,0.2797053915,0.1294849662);
  EPS=1.0e-12;
  K=7;
var
  S,P,H,AA,BB,W,X,Q,TempF1,TempF2:extended;
  G:T3DVector;  
  I,J,M,N,L:integer;
begin
  M:=30;  //can be changed
  for L:=1 to 3 do
  begin
    TempF1:=CmplxABS(F1[L]);
    TempF2:=CmplxABS(F2[L]);
    S:=(B-A)*0.001;
    P:=0.0;
    repeat
      H:=(B-A)/M;
      G[L]:=0.0;
      for I:=1 to M do
      begin
        AA:=A+(I-1)*H;
        BB:=A+I*H;
        W:=0.0;
        for J:=1 to K do
        begin
          X:=((BB-AA)*IntegPoint[J]+(BB+AA))/2.0;
          W:=W+sqr(TempF1*sin(2*Pi*Frequency[0]*X-Pi/4.0)+TempF2*sin(2*Pi*Frequency[1]*X-Pi/4.0))*coef[J];
        end;
        G[L]:=G[L]+W;
      end;
      G[L]:=G[L]*H/2.0;
      Q:=abs(G[L]-P)/(1.0+abs(G[L]));
      P:=G[L];
      Inc(M);
    until (Q<EPS)or(abs(H)<=abs(S));
  end;
  Result:=Model(G);
end;

end.