/*-----------------------------------------------------------------------*/
/* RC circular section (N,M,S) */
/*-----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define pi M_PI
void main_part(char fnameR[],char fnameW[]);
double fnc(double alpha,double n,double p,double e,double r,double rs);
/*-----------------------------------------------------------------------*/

int main(int argc,char *argv[])
{
    char fnameR[50]; /*入力ファイル名*/
    char fnameW[50]; /*出力ファイル名*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);

    main_part(fnameR,fnameW);

    printf("リターンキーを押してください");
    getchar();
    return 0;
}
/*---------------------------------------------------------------------------*/
void main_part(char fnameR[],char fnameW[])
{
    FILE *fin,*fout;
    char dat[256];
    char sv[20];
    double r,rs,As,n;
    double NN,MM,SS;
    double AI,II,aa,theta;
    double e,p;
    double a1,a2,ai,f1,f2,fi,alpha,yn;
    double q1;
    double sig_cc,sig_cd,sig_ss,sig_sd,tau_c;

    /* データ読み込み */
    fin=fopen(fnameR,"r");
    fout=fopen(fnameW,"w");
    fprintf(fout,"RC circular section (N,M,S)\n");
    fprintf(fout,"Input data,,,,,,,Compressive side,,Tensile side,,Average,Neutral axis\n");
    fprintf(fout,"r,rs,As,n,N,M,S,sig_c,sig_s',sig_c',sig_s,tau_c,alpha(deg.),yn\n");
    while(fgets(dat,sizeof dat,fin)!=NULL){
        r =atof(strtok(dat,","));  /*断面半径*/
        rs=atof(strtok(NULL,",")); /*鉄筋位置半径*/
        As=atof(strtok(NULL,",")); /*鉄筋断面積*/
        n =atof(strtok(NULL,",")); /*弾性係数比*/
        NN=atof(strtok(NULL,",")); /*軸力*/
        MM=atof(strtok(NULL,",")); /*モーメント*/
        SS=atof(strtok(NULL,",")); /*せん断力*/

        alpha=0.0;yn=0.0;
        sig_cc=0.0;sig_cd=0.0;
        sig_ss=0.0;sig_sd=0.0;
        AI=pi*r*r+n*As;
        II=pi*r*r*r*r/4.0+n*As*rs*rs/2.0;
        sig_cc=NN/AI+MM/II*r; /*圧縮縁コンクリート応力*/
        sig_cd=NN/AI-MM/II*r; /*引張縁コンクリート応力*/
        sig_ss=rs/r*n*sig_cd; /*引張縁鉄筋*/
        sig_sd=rs/r*n*sig_cc; /*圧縮縁鉄筋応力*/

        sig_cc=-sig_cc;
        sig_cd=-sig_cd;
        sig_ss=-sig_ss;
        sig_sd=-sig_sd;

        if(0.0<sig_cd){
            alpha=0.0;yn=0.0;
            sig_cc=0.0;sig_cd=0.0;
            sig_ss=0.0;sig_sd=0.0;
            e=MM/NN;
            p=As/pi/r/r;
            a1=0.0;
            a2=pi;
            ai=(a1+a2)/2.0;
            do{
                f1=fnc(a1,n,p,e,r,rs);
                f2=fnc(a2,n,p,e,r,rs);
                fi=fnc(ai,n,p,e,r,rs);
                if(fi==0.0)break;
                if(f1*fi<0.0)a2=ai;
                if(fi*f2<0.0)a1=ai;
                ai=(a1+a2)/2.0;
            }while(1e-6<fabs(fi));
            alpha=ai;
            q1=2.0/3.0*sin(alpha)*sin(alpha)*sin(alpha)+sin(alpha)*cos(alpha)*cos(alpha)-alpha*cos(alpha);
            sig_cc=NN*(1.0-cos(alpha))/r/r/(q1-n*pi*p*cos(alpha));
            sig_ss=(rs/r+cos(alpha))/(1.0-cos(alpha))*n*sig_cc;
            sig_sd=(rs/r-cos(alpha))/(1.0-cos(alpha))*n*sig_cc;

            sig_cc=-sig_cc;
            sig_sd=-sig_sd;

            yn=r*cos(pi-alpha);
        }
        theta=acos(rs/r);
        aa=pi*r*r-(r*r*theta-rs*r*sin(theta));
        tau_c=SS/aa;

              printf("%.1f %.1f %.1f %.0f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.4f %.3f %.3f\n",r,rs,As,n,NN,MM,SS,sig_cc,sig_sd,sig_cd,sig_ss,tau_c,alpha/pi*180.0,yn);
        fprintf(fout,"%.1f,%.1f,%.1f,%.0f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.4f,%.3f,%.3f\n",r,rs,As,n,NN,MM,SS,sig_cc,sig_sd,sig_cd,sig_ss,tau_c,alpha/pi*180.0,yn);
    }
    fclose(fout);
    fclose(fin);


}
/*---------------------------------------------------------------------------*/
double fnc(double alpha,double n,double p,double e,double r,double rs)
{
    double q1,q2,ff;

    q1=2.0/3.0*sin(alpha)*sin(alpha)*sin(alpha)+sin(alpha)*cos(alpha)*cos(alpha)-alpha*cos(alpha);
    q2=alpha/4.0-1.0/6.0*sin(alpha)*sin(alpha)*sin(alpha)*cos(alpha)-1.0/8.0*sin(2.0*alpha);
    ff=(q1-n*pi*p*cos(alpha))*e/r-(q2+n*pi*p/2.0*(rs*rs/r/r));

    return ff;
}
/*---------------------------------------------------------------------------*/