/*---------------------------------------------------------------------------*/
/* Amstutzによる座屈圧力計算式 */
/*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#define pi M_PI

void calc1(double Ess,double rm,double t,double SigN,double *Phi,double *Psi,double *Omega,double *theta);
double calc2(double Ess,double rm,double t,double k0,double m,double SigF,double SigN,double Phi,double Psi);
/*---------------------------------------------------------------------------*/
int main()
{
    FILE *fp;
    double Es,Eps,Pos,As,Bg;
    double D0,t0,Sigy,Siga,Jc,Temp;
    double Ess,ny,SigF,t,rm,r0d,k0,SigN,Pk,offset,DeltaD,m;
    double Phi,Psi,Omega,theta;
    double c1,c2,c3,x1,x2,x3;
    double pkr[5];

    /* 座屈荷重計算(補剛材無し) */
    /* D0:鉄管内径 */
    /* t0:設計板厚 */
    /* Eps:鋼材の余裕厚 */
    /* Sigy:鋼材の降伏点 */
    /* Siga:鋼材の許容応力 */
    /* Pos:鋼材のﾎﾟｱｿﾝ比 */
    /* Temp:鉄管の温度降下量 */
    /* Jc:継手効率 */
    /* Pk:座屈圧力 */
    /* offset:段違い */
    /* DeltaD:初期不整 */

    int nnn,iii;

    Es=206000;
    Eps=1.5;
    Pos=0.3;
    As=1.2e-5;
    Temp=20.0;
    Jc=1.00;
    Bg=1.0;
    t0=30.0;

    fp=fopen("inp_PS_AMSid.prn","w");
    fprintf(fp,"#D0  t0  rm/t0  Pk1  Pk2  Pk3  Pk4  Pk5\n");
    /*==========================================*/
    for(nnn=35;nnn<=140;nnn++){
    /*==========================================*/

        D0=2.0*(double)nnn*t0;

        for(iii=0;iii<=4;iii++){
            switch(iii){
                case 0:Sigy=885.0;Siga=400.0;break;/*HT100*/
                case 1:Sigy=685.0;Siga=330.0;break;/*SHY685*/
                case 2:Sigy=450.0;Siga=240.0;break;/*SM570*/
                case 3:Sigy=315.0;Siga=175.0;break;/*SM490B*/
                case 4:Sigy=235.0;Siga=130.0;break;/*SM400B*/
            }

            DeltaD=0.0;
            offset=0.0;

            DeltaD=D0/300.0;
            offset=0.05*t0;

            Ess=Es/(1.0-Pos*Pos);
            ny=1.5-0.5*1.0/(1.0+0.002*Es/Sigy)/(1.0+0.002*Es/Sigy);
            SigF=ny*Sigy/sqrt(1.0-Pos+Pos*Pos);

            t=t0-Eps;
            rm=(D0+t0)/2.0;
            r0d=(D0+2.0*t0)/2.0;
            k0=(As*Temp+Bg*Siga*Jc/Es)*r0d/(1.0+Bg);
            /*k0=0.4e-3*rm;

            /*初期不整の考慮*/
            rm=rm*(1.0+1.52240775*DeltaD/D0);
            r0d=r0d*(1.0+1.52240775*DeltaD/D0);

            /*段違い(offset)の考慮*/
            m=1.0+3.0*offset/t;

            /*二分法によるSigNの計算*/
            x1=Ess*(3.0*3.0-1.0)/12.0*t*t/rm/rm;
            x2=SigF;
            do{
                x3=0.5*(x1+x2);
                SigN=x1;calc1(Ess,rm,t,SigN,&Phi,&Psi,&Omega,&theta);c1=calc2(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi);
                SigN=x2;calc1(Ess,rm,t,SigN,&Phi,&Psi,&Omega,&theta);c2=calc2(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi);
                SigN=x3;calc1(Ess,rm,t,SigN,&Phi,&Psi,&Omega,&theta);c3=calc2(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi);
                if(c1*c3<=0.0){x1=x1;x2=x3;}
                if(c3*c2<0.0) {x1=x3;x2=x2;}
                if(c1*c2>0.0){
                    printf("Imposible  x1=%.3f  x2=%.3f  c1=%.3f  c2=%.3f\n",x1,x2,c1,c2);
                    return 0;
                }
            }while(fabs(x1-x2)>1e-6);

            Pk=SigN/(rm/t)/(1.0+2.0*Omega*rm/t*(SigF-m*SigN)/Ess);
            pkr[iii]=Pk;

            /*printf("--------------------------------------------\n");*/
            /*printf("D0=%.1f t0=%.1f Sigy=%.1f Siga=%.1f k0=%.3f offset=%.3f DeltaD=%.3f\n",D0,t0,Sigy,Siga,k0,offset,DeltaD);*/
            /*printf("SigN=%.3f SigF=%.3f Pk=%.3f\n",SigN,SigF,Pk);*/
            /*printf("Phi=%.3f Psi=%.3f Omega=%.3f ang=%.3f\n",Phi,Psi,Omega,theta);*/
        }
        printf("%.1f %.1f %.1f %.3f %.3f %.3f %.3f %.3f\n",D0,t0,0.5*D0/t0,pkr[0],pkr[1],pkr[2],pkr[3],pkr[4]);
        fprintf(fp,"%.1f %.1f %.1f %.3f %.3f %.3f %.3f %.3f\n",D0,t0,0.5*D0/t0,pkr[0],pkr[1],pkr[2],pkr[3],pkr[4]);
    /*==========================================*/
    }
    /*==========================================*/
    fclose(fp);

    return 0;
}
/*---------------------------------------------------------------------------*/
void calc1(double Ess,double rm,double t,double SigN,double *Phi,double *Psi,double *Omega,double *theta)
{
    double ang,cc,ee,alpha,delta,beta,cota,gamma;
    /* ang=90〜0度 */
    /* ee*ang=270〜180度 */

    ee=sqrt(1.0+12.0*rm*rm/t/t*SigN/Ess);
    ang=270.0/ee;
    do{
        ang=ang-0.0001;
        alpha=ang*pi/180.0;
        cc=ee*tan(alpha)-tan(ee*alpha);
    }while(cc<0.0);
    delta=(ee*ee-1.0)*(1.0-cos(ee*alpha));
    beta=(ee-1.0/ee)*(ee*alpha*cos(ee*alpha)-sin(ee*alpha));
    if(ang==90.0){cota=0.0;}else{cota=1.0/tan(alpha);}
    gamma=-ee*(ee*alpha-sin(ee*alpha)*cos(ee*alpha)+ee*alpha*sin(ee*alpha)*sin(ee*alpha)/sin(alpha)/sin(alpha)-ee*sin(ee*alpha)*sin(ee*alpha)*cota);

    *Phi=ee*ee*ee*beta/pi/delta;
    *Psi=-gamma/4.0/beta/delta;
    *Omega=-cos(ee*alpha)/(1.0-cos(ee*alpha));
    *theta=alpha/pi*180.0;
}
/*---------------------------------------------------------------------------*/
double calc2(double Ess,double rm,double t,double k0,double m,double SigF,double SigN,double Phi,double Psi)
{
    double a,b;

    a=(k0/rm+SigN/Ess)*pow(1.0+12.0*rm*rm/t/t*SigN/Ess,1.5);
    b=2.0*Phi*rm/t*(SigF-m*SigN)/Ess*(1.0-2.0*Psi*rm/t*(SigF-m*SigN)/Ess);

return a-b;
}
/*---------------------------------------------------------------------------*/