11楼: Originally posted by dingd at 2012-03-26 15:14:45:
代码中有“M/abs(M)”，而M有被赋值为0的时候（注意yy数组中有0项），这样会导致被0除的错误。

Title "Type your title here"
Constant xx=[0.703,1.698,1.698,0.703,-0.703,-1.698,-1.698,-0.703,1.246,3.008,3.008,1.246,-1.246,
             -3.008,-3.008,-1.246,0,-2.926,-4.138,-2.926,0,2.926,4.138,2.926],
         yy=[1.698,0.703,-0.703,-1.698,-1.698,-0.703,0.703,1.698,3.008,1.246,-1.246,-3.008,-3.008,
            -1.246,1.246,3.008,4.138,2.926,0,-2.926,-4.138,-2.926,0,2.926];
Parameter a[1,4], b[1,4],theta[0,pi/2],phi[-pi/8,pi/8];
Minimum;
StartProgram;
var i,j: Integer;
    Etot,r,angle,out1,out2,out3: double;
    Ax,Ay,Ar,Ae,Bx,By,Br,Be,dist,rou: double;
    M: array[0..3] of array[1..24] of array[1..4] of double;
Begin
        for i:=1 to 24 do begin
            M[0,i,1]:=xx[i];
            M[0,i,2]:=yy[i];
        end;
        for i:=1 to 8 do begin
            M[0,i,3]:= 1.9;
            M[0,i,4]:= 0.044;
            M[0,i+8,3]:= 1.9;
            M[0,i+8,4]:= 0.044;
            M[0,i+16,3]:= 2.11;
            M[0,i+16,4]:= 0.202;
        end;
        for i := 1 to 24 do
            for j := 1 to 4 do begin
                M[1,i,j]:=M[0,i,j];
                M[2,i,j]:=M[0,i,j];
                M[3,i,j]:=M[0,i,j];
            end;
        for i:=1 to 24 do begin
            r:=sqrt(M[1,i,1]*M[1,i,1]+M[1,i,2]*M[1,i,2]);
            if M[1,i,1]>=0 then
               angle:=arcsin(M[1,i,2]/r)
            else begin
               //angle:= (M[1,i,2]/abs(M[1,i,2]))*(pi-abs(arcsin(M[1,i,2]/r)));
               if M[1,i,2] >=0 then
                  angle:= (pi-abs(arcsin(M[1,i,2]/r)))
               else
                  angle:= -(pi-abs(arcsin(M[1,i,2]/r)));
            end;

            angle:=angle+phi;
            M[1,i,1]:=r*cos(angle);
            M[1,i,2]:=r*sin(angle);
        end;

        for i:=1 to 24 do begin
            M[2,i,1]:=M[1,i,1]+b;
            M[3,i,1]:=M[1,i,1]+a*cos(theta);
            M[3,i,2]:=M[1,i,2]+a*sin(theta);
        end;
        Etot:=0;
        for i:=1 to 24 do
            for j:=1 to 24 do begin
                Ax:=M[1,i,1];
                Ay:=M[1,i,2];
                Ar:=M[1,i,3];
                Ae:=M[1,i,4];
                Bx:=M[2,j,1];
                By:=M[2,j,2];
                Br:=M[2,j,3];
                Be:=M[2,j,4];
                dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By));
                rou:=dist/(Ar+Br);
                if dist>3.311 then
                   out1:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou)))
                else
                   out1:= 336.176*sqrt(Ae*Be)/(rou*rou);
                Bx:=M[3,j,1];
                By:=M[3,j,2];
                Br:=M[3,j,3];
                Be:=M[3,j,4];
                dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By));
                rou:=dist/(Ar+Br);
                if dist>3.311 then
                   out2:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou)))
                else
                   out2:= 336.176*sqrt(Ae*Be)/(rou*rou);
                Ax:=M[2,i,1];
                Ay:=M[2,i,2];
                Ar:=M[2,i,3];
                Ae:=M[2,i,4];
                dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By));
                rou:=dist/(Ar+Br);
                if dist>3.311 then
                   out3:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou)))
                else
                   out3:= 336.176*sqrt(Ae*Be)/(rou*rou);
                Etot:=out1+out2+out3+Etot;
            end;
        FunctionResult:= Etot;
End;
EndProgram;

14楼: Originally posted by dingd at 2012-03-26 15:41:42:
1stOpt代码：

Title "Type your title here"
Constant xx=,
         yy=;
Parameter a, b,theta,phi;
Minimum;
StartProgram;
var i,j: Integer;
    Etot,r,angle,out1,out2,out3: double; ...