瑞典条分法计算代码(恳请指正)

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<malloc.h>

void A_Parameter(void);
void Coordinate(void); //该函数主要用于生成圆心坐标
void Arc(double R,double AngleTime,double DotX);
void GravityComponet(double *NetSoilStripWeight,double *ArcLength,double *ArcAngle,double r,double AngleTime);
void SafetyFactor(double MR,double MT,double R,double OneAngleTime);

#define Angle 180/(4*atan(1))
#define radian 4*atan(1)/180

double H_Distance=0;
int SoilStripNumber=0;     //土条数
int S_L_Number=0;          //将土层层数定义为全局变量
int numberone=0;

struct result{
    double Rdius;
    double Ks;
    double AngleTime;
    double Mr;    //抗滑力矩；
    double Mt;    //下滑力矩；
}Result;

struct AverageParameter{
     double VolumWeight;
     double FrictionAngle;
     double Cohesive;
     double Height;
}Average;

int main()
{ 
    A_Parameter();
    printf("输出加权平均后的土壤参数(容重 内摩擦角 黏聚力)和坡高\n");
    printf("%lf  %lf  %lf  %lf\n", Average .VolumWeight, Average .FrictionAngle, Average .Cohesive, Average.Height);
    Coordinate();       //该函数主要用于生成圆心坐标
    //printf("numberone is %d\n",numberone);0
    printf("输出结果\n");
    printf("根据瑞典条分法求得边坡安全系数K:  K=%.2lf\n",Result.Ks);
    printf("根据瑞典条分法求得滑移面半径R:  R=%.2lf\n",Result.Rdius);
printf("根据瑞典条分法求得圆心坐标(X,Y):  X=%.2lf  Y=%.2lf\n",-Result.Rdius*cos(Result.AngleTime*radian),Result.Rdius
*sin(Result.AngleTime*radian));
    printf("根据瑞典条分法求得抗滑力矩Mr:  Mr=%.2lf\n",Result.Mr);
    printf("根据瑞典条分法求得下滑力矩Mt:  Mt=%.2lf\n",Result.Mt);
    printf("根据瑞典条分法求得半径与水平面的夹角Angle:  Angle=%.2lf\n",Result.AngleTime);
   
    return 0;
}
void A_Parameter ()
{
    double *VolumWeight,*FrictionAngle;
    double *Cohesive,*Thickness;
int i=0,j=0;

    printf("输入要划分的土条数\n");
    scanf("%d",&SoilStripNumber);
    printf("输入边坡水平距离\n");
    scanf("%lf",&H_Distance);
    printf("输入土层层数\n");
scanf("%d",&S_L_Number);

    VolumWeight=(double *)calloc(S_L_Number,sizeof(double));
    FrictionAngle=(double *)calloc(S_L_Number,sizeof(double));
    Cohesive=(double *)calloc(S_L_Number,sizeof(double));
Thickness=(double *)calloc(S_L_Number,sizeof(double));

     printf("输入各土层容重\n");
    while(i<S_L_Number){
        scanf("%lf",&(VolumWeight[i]));
        i++;
    }
    
    printf("输入各土层内摩擦角单位为度\n");
    for(j=0;j<S_L_Number;j++)
        scanf("%lf",&(FrictionAngle[j]));
       
    printf("输入各土层的黏聚力单位为千帕\n");   
    for(j=0;j<S_L_Number;j++)
        scanf("%lf",&(Cohesive[j]));
        
    printf("输入各土层的厚度单位为米\n");
    for(j=0;j<S_L_Number;j++)
        scanf("%lf",&(Thickness[j]));
        
    double T_VolumWeight=0,T_FrictionAngle=0;  //各层高度和各层参数乘积之和。
    double T_Cohesive=0,T_Thickness=0;         //各黏聚力和各层高度乘积之和。总的厚度。
    
    for(j=0;j<S_L_Number;j++)
        T_VolumWeight+=VolumWeight[j]*Thickness[j];      //计算各层重度与各层高度乘积之和
    
    for(j=0;j<S_L_Number;j++)
        T_FrictionAngle+=FrictionAngle[j]*Thickness[j];  //计算各层内摩擦角和各层高度乘积之和
        
    for(j=0;j<S_L_Number;j++)
        T_Cohesive+=Cohesive[j]*Thickness[j];  //计算各层黏聚力与各层高度乘积之和
   
    for(j=0;j<S_L_Number;j++)
        T_Thickness+=Thickness[j];  //计算总得土层厚度
    
    Average .VolumWeight=T_VolumWeight/T_Thickness;
    Average .FrictionAngle=T_FrictionAngle/T_Thickness;
    Average .Cohesive=T_Cohesive/T_Thickness;
    Average .Height=T_Thickness;
    
    free(VolumWeight);
    free(FrictionAngle);
    free(Cohesive);
    free(Thickness);
}
void Coordinate(void)
{
int i=0,j=0;
double a;
    double DotX=0,DotX1=0,AngleOne=0;;
    double DeltaR=0,R=0,AngleTime=0;
    DeltaR=(3*Average.Height)/200;
  
    DotX1=H_Distance;
    for(i=0;i<=200 ;i++){
        R=Average.Height+i*DeltaR;
        AngleOne=(asin(sqrt(pow(H_Distance,2)+pow(Average.Height,2))/(2*R))+atan(DotX1/Average.Height))*Angle;
        a=(90-AngleOne)/ SoilStripNumber;  
        //printf("\n");
        for(j=0;j<SoilStripNumber;j++){
            AngleTime=(AngleOne+j*a);
            DotX=sqrt(pow(R,2)-pow((Average.Height-(R*sin(AngleTime*radian))),2))-R*cos(AngleTime*radian);
            if((R*sin(AngleTime*radian)>=Average.Height)&&(DotX<=Average.Height)){
                //printf("  %.2lf",DotX);
                Arc (R,AngleTime,DotX);
            }
        }
    }

}
void Arc(double R,double AngleTime,double DotX)  
{
    int i=0;
    double AngleTimeOne=0,DeltaAngleTimeOne=0,AngleTimetwo=0;
    double Distance=0;
    double r=R,OneAngleTime=AngleTime;     //OneAngleTime做函数的实参；
    double *S_LineFy=0;  //直线部分的竖距
    double *ArcLength,*ArcAngle,*ArcX;
    double *ArcFy,*ArcSoilStripWeight;
    double *S_LineSoilStripWeight;
    double *NetSoilStripWeight;
   
    S_LineFy=(double *)calloc(SoilStripNumber,sizeof(double));
    ArcLength=(double *)calloc(SoilStripNumber,sizeof(double));
    ArcAngle=(double *)calloc(SoilStripNumber,sizeof(double));
    ArcX=(double *)calloc(SoilStripNumber,sizeof(double));
    ArcFy=(double *)calloc(SoilStripNumber,sizeof(double));
    ArcSoilStripWeight=(double *)calloc(SoilStripNumber,sizeof(double));
    S_LineSoilStripWeight=(double *)calloc(SoilStripNumber,sizeof(double));
    NetSoilStripWeight=(double *)calloc(SoilStripNumber,sizeof(double));
    
    if((r*sin(AngleTime*radian))>=Average.Height){
        Distance=r*cos(AngleTime*radian)+DotX;   //Distance表示水平距离
        AngleTimetwo=acos(Distance/r)*Angle;
        DeltaAngleTimeOne=(AngleTime-AngleTimetwo)/SoilStripNumber;
    }
    
    //printf("  %.2lf",DotX);
    //竖条高和竖条自重
    for(i=0;i<SoilStripNumber;i++){
        AngleTimeOne=DeltaAngleTimeOne*(i+1);
        //printf("  %.2lf",SoilStripWide);
        if(AngleTimeOne<=AngleTime){
            ArcFy[i]=r*sin(AngleTime*radian)-r*cos(((90-AngleTime)+AngleTimeOne)*radian);
            ArcX[i]=r*sin(((90-AngleTime)+AngleTimeOne)*radian)-r*cos(AngleTime*radian);
            //printf("  %.2lf",ArcFy[i]);
        }
        //printf("  %.2lf",ArcX[i]);
        if(i==1){
            ArcSoilStripWeight[i]=Average.VolumWeight*ArcX[i]*ArcFy[i];
        }else{
            ArcSoilStripWeight[i]=Average.VolumWeight*(ArcX[i]-ArcX[i-1])*ArcFy[i];
        }
        //printf("  %.2lf",ArcSoilStripWeight[i]);
    }
    //printf("\n");
    for(i=0;i<SoilStripNumber;i++){
        
        if(ArcX[i]<=H_Distance){
            if(i==0){
                S_LineFy[i]=(Average.Height/H_Distance)*ArcX[i];
                S_LineSoilStripWeight[i]=Average.VolumWeight*ArcX[i]*S_LineFy[i];
            }
            else{
                S_LineFy[i]=(Average.Height/H_Distance)*ArcX[i];
                S_LineSoilStripWeight[i]=Average.VolumWeight*(ArcX[i]-ArcX[i-1])*S_LineFy[i];
            }
        }
        else{
            S_LineFy[i]=Average.Height;
            S_LineSoilStripWeight[i]=Average.VolumWeight*(ArcX[i]-ArcX[i-1])*Average.Height;
        }
    //printf("  %.2lf",S_LineSoilStripWeight[i]);
    }
    
    //计算直线与圆弧所对应的土条重量差值,弧长,角度
    for(i=0;i<SoilStripNumber;i++){
        if(i<(SoilStripNumber-1))
            NetSoilStripWeight[i]=S_LineSoilStripWeight[i]-ArcSoilStripWeight[i];
        else
            NetSoilStripWeight[i]=0.5*(S_LineSoilStripWeight[i-1]-ArcSoilStripWeight[i-1]);
    //printf("  %.2lf",NetSoilStripWeight[i]);
    }
    //printf("\n");   
    for(i=0;i<SoilStripNumber;i++){
        if(i==0){
            ArcLength[i]=sqrt(pow(ArcX[i],2)+pow(ArcFy[i],2));
            ArcAngle[i]=atan2(ArcFy[i],ArcX[i]);  //弧度制
        }
        else{
            ArcLength[i]=sqrt(pow(ArcX[i]-ArcX[i-1],2)+pow((ArcFy[i]-ArcFy[i-1]),2));
            ArcAngle[i]=atan2((ArcFy[i]-ArcFy[i-1]),(ArcX[i]-ArcX[i-1]));  //弧度制
        }
    //printf("  %.2lf",ArcAngle[i]);
    }
    
    GravityComponet(NetSoilStripWeight,ArcLength,ArcAngle,r,OneAngleTime);
    free(S_LineFy);
    free(ArcLength);
    free(ArcAngle);
    free(ArcX);
    free(ArcFy);
    free(ArcSoilStripWeight);
    free(S_LineSoilStripWeight);
    free(NetSoilStripWeight);
}

void GravityComponet(double *NetSoilStripWeight,double *ArcLength,double *ArcAngle,double r,double AngleTime)  //计算下滑力和抗滑力
{
    int i=0;
    double MR=0,MT=0,R=r,OneAngleTime=AngleTime;  //抗滑力矩MT和滑动力矩MR
    double *TangantGravity,*VerticalGravity;
    double *ShearStrength,*ShearForce;            //微段弧长上的剪应力和剪力
    
    TangantGravity=(double *)calloc(SoilStripNumber,sizeof(double));
    VerticalGravity=(double *)calloc(SoilStripNumber,sizeof(double));
    ShearStrength=(double *)calloc(SoilStripNumber,sizeof(double));
    ShearForce=(double *)calloc(SoilStripNumber,sizeof(double));
    //计算下滑力
    for(i=0;i<SoilStripNumber;i++){
        TangantGravity[i]=NetSoilStripWeight[i]*sin(ArcAngle[i]);
        VerticalGravity[i]=NetSoilStripWeight[i]*cos(ArcAngle[i]);     //垂直于滑动面的法向力
    }
    //计算抗剪强度
    for(i=0;i<SoilStripNumber;i++){
        ShearStrength[i]=(VerticalGravity[i]/ArcLength[i])*tan(Average.FrictionAngle*radian)+Average.Cohesive;
        ShearForce[i]=ShearStrength[i]*ArcLength[i];
    }
    //printf("\n");
    for(i=0;i<SoilStripNumber;i++){
        MT+=(TangantGravity[i]*r);
        MR+=(ShearForce[i]*r);
    }
    
    //printf("  %.2lf",MT);
    SafetyFactor(MR,MT,R,OneAngleTime);     //计算安全系数
    free(TangantGravity);
    free(VerticalGravity);
    free(ShearStrength);
    free(ShearForce);
}
void SafetyFactor(double MR,double MT,double R,double OneAngleTime)
{
    numberone++;
    double Ks_one=0;
    Ks_one=MR/MT;
    if(numberone==1){
        Result.Ks=Ks_one;
        Result.Rdius=R;
        Result.AngleTime=OneAngleTime;
        Result.Mr=MR;
        Result.Mt=MT;
    }   
    //下面代码将最小的Ks_one赋值给Ks
    if(numberone!=1){
        if(Result.Ks>=Ks_one){
            Result.Ks=Ks_one;
            Result.Rdius=R;
            Result.AngleTime=OneAngleTime;
            Result.Mr=MR;
            Result.Mt=MT;
        }
    }
}