#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 99
#define Mmax 4

void main_part(char fnameR[],char fnameW[]);
void STIF(int STRES,double EE,double po,double *AA,double *BB);
void CALC0(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb);
void CALC1(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb);
void CALC2(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb);
void CALC3(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb);
void MATGJ(int n,double a[Mmax+1][Mmax+2]);
/*----------------------------------------------------------------------------*/

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];
    
    int k,iii,jjj;
    char strcom[256]; /* コメント文 */
    int nd;           /* 計算組数 */
    int STRES;        /* 0;平面歪み，1：平面応力 */
    double P[Nmax+1];        /* 内圧 */
    double a[Nmax+1];        /* 水路内半径 */
    double b[Nmax+1];        /* 掘削半径 */
    double Ec[Nmax+1];       /* コンクリート弾性係数 */
    double nc[Nmax+1];       /* コンクリートポアソン比 */
    double Es[Nmax+1];       /* 鉄筋弾性係数 */
    double ta[Nmax+1];       /* 内側鉄筋等価板厚 */
    double tb[Nmax+1];       /* 外側鉄筋等価板厚 */
    double Eg[Nmax+1];       /* 岩盤弾性係数 */
    double ng[Nmax+1];       /* 岩盤ポアソン比 */
    double sigsa[Nmax+1][4]; /* 内側鉄筋応力 */
    double sigsb[Nmax+1][4]; /* 外側鉄筋応力 */
    double sigrc[Nmax+1][4]; /* コンクリート半径方向応力 */
    double sigtc[Nmax+1][4]; /* コンクリート円周方向応力 */
    double ua[Nmax+1][4];    /* 内壁半径方向変位 */
    double ub[Nmax+1][4];    /* 外壁半径方向変位 */

    /*作業用変数 */
    double PP,ar,br;
    double Ecc,ncc;
    double Ess,taa,tbb;
    double Egg,ngg;
    double ssa,ssb;
    double src,stc;
    double uaa,ubb;

    /* テキスト入力 */
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);strcpy(strcom,strtok(dat,","));
    fgets(dat,sizeof dat,fin);nd=atoi(strtok(dat,","))-1;STRES=atoi(strtok(NULL,","));
    k=0;
    while(fgets(dat,sizeof dat,fin)!=NULL){
        P[k] =atof(strtok(dat,","));
        a[k] =atof(strtok(NULL,","));
        b[k] =atof(strtok(NULL,","));
        Ec[k]=atof(strtok(NULL,","));
        nc[k]=atof(strtok(NULL,","));
        Es[k]=atof(strtok(NULL,","));
        ta[k]=atof(strtok(NULL,","));
        tb[k]=atof(strtok(NULL,","));
        Eg[k]=atof(strtok(NULL,","));
        ng[k]=atof(strtok(NULL,","));
        k=k+1;
    }
    fclose(fin);
    nd=k-1;

    for(iii=0;iii<=nd;iii++){
        PP=P[iii];ar=a[iii];br=b[iii];
        Ecc=Ec[iii];ncc=nc[iii];
        Ess=Es[iii];taa=ta[iii];tbb=tb[iii];
        Egg=Eg[iii];ngg=ng[iii];
        for(jjj=0;jjj<=3;jjj++){
            switch(jjj){
                case 0: /* ひび割れのない無筋コンクリート */
                    CALC0(STRES,PP,ar,br,Ecc,ncc,Ess,taa,tbb,Egg,ngg,&ssa,&ssb,&src,&stc,&uaa,&ubb);break;
                case 1: /* ひび割れの無い単鉄筋コンクリート */
                    CALC1(STRES,PP,ar,br,Ecc,ncc,Ess,taa,tbb,Egg,ngg,&ssa,&ssb,&src,&stc,&uaa,&ubb);break;
                case 2: /* ひび割れのある単鉄筋コンクリート */
                    CALC2(STRES,PP,ar,br,Ecc,ncc,Ess,taa,tbb,Egg,ngg,&ssa,&ssb,&src,&stc,&uaa,&ubb);break;
                case 3: /* ひび割れのある複鉄筋コンクリート */
                    CALC3(STRES,PP,ar,br,Ecc,ncc,Ess,taa,tbb,Egg,ngg,&ssa,&ssb,&src,&stc,&uaa,&ubb);break;
            }
            sigsa[iii][jjj]=ssa; /* 内側鉄筋応力 */
            sigsb[iii][jjj]=ssb; /* 外側鉄筋応力 */
            sigrc[iii][jjj]=src; /* コンクリート半径方向応力 */
            sigtc[iii][jjj]=stc; /* コンクリート円周方向応力 */
            ua[iii][jjj]=uaa;    /* 内壁半径方向変位 */
            ub[iii][jjj]=ubb;    /* 外壁半径方向変位 */
        }
    }

    /* テキスト出力 */
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s\n",strcom);
    fprintf(fout,"nd,%d,STRES,%d\n",nd+1,STRES);
    fprintf(fout,"入力データ\n");
    fprintf(fout,"i,P,a,b,Ec,nc,Es,ta,tb,Eg,ng\n");
    for(iii=0;iii<=nd;iii++){
        sprintf(sv,"%d",iii+1)   ;strcpy(dat,sv);
        sprintf(sv,"%e", P[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e", a[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e", b[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Ec[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",nc[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Es[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ta[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",tb[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Eg[iii]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ng[iii]) ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    for(jjj=0;jjj<=3;jjj++){
        switch(jjj){
            case 0:strcpy(dat,"ひび割れのない無筋コンクリート")  ;break;
            case 1:strcpy(dat,"ひび割れのない単鉄筋コンクリート");break;
            case 2:strcpy(dat,"ひび割れのある単鉄筋コンクリート");break;
            case 3:strcpy(dat,"ひび割れのある複鉄筋コンクリート");break;
        }
        fprintf(fout,"%s\n",dat);
        fprintf(fout,"i,σsa,σsb,σrc,σtc,u(r=a),u(r=b)\n");
        for(iii=0;iii<=nd;iii++){
            sprintf(sv,"%d",iii+1)           ;strcpy(dat,sv);
            sprintf(sv,"%e",sigsa[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%e",sigsb[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%e",sigrc[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%e",sigtc[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%e",   ua[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            sprintf(sv,"%e",   ub[iii][jjj]) ;strcat(dat,",");strcat(dat,sv);
            fprintf(fout,"%s\n",dat);
        }
    }
    fclose(fout);
}
/*----------------------------------------------------------------------------*/

void STIF(int STRES,double EE,double po,double *AA,double *BB)
{
    switch(STRES){
        case 0: /* 平面ひずみ */
            *AA=EE*(1.0-po)/(1.0+po)/(1.0-2.0*po);
            *BB=EE*po/(1.0+po)/(1.0-2.0*po);
            break;
        case 1: /* 平面応力 */
            *AA=EE/(1.0-po*po);
            *BB=EE*po/(1.0-po*po);
            break;
    }
}
/*----------------------------------------------------------------------------*/

void CALC0(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb)
{
    /* ひび割れのない無筋コンクリート */
    int mm=2;
    double sm[Mmax+1][Mmax+2];
    double Ac,Bc,C1,C2,Cg;
    int i,j;

    STIF(STRES,Ecc,ncc,&Ac,&Bc);

    for(i=0;i<=mm;i++){
        for(j=0;j<=mm+1;j++){
            sm[i][j]=0.0;
    }}
    sm[0][0]=Ac+Bc;
    sm[0][1]=-(Ac-Bc)/ar/ar;
    sm[1][0]=br;
    sm[1][1]=1.0/br;
    sm[1][2]=-1.0/br;
    sm[2][0]=Ac+Bc;
    sm[2][1]=-(Ac-Bc)/br/br;
    sm[2][2]=Egg/(1.0+ngg)/br/br;
    sm[0][3]=-PP;
    MATGJ(mm,sm);
    C1=sm[0][3];
    C2=sm[1][3];
    Cg=sm[2][3];
    *ssa=0.0;
    *ssb=0.0;
    *src=(Ac+Bc)*C1-(Ac-Bc)*C2/ar/ar;
    *stc=(Ac+Bc)*C1+(Ac-Bc)*C2/ar/ar;
    *uaa=C1*ar+C2/ar;
    *ubb=C1*br+C2/br;

}
/*----------------------------------------------------------------------------*/

void CALC1(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb)
{
    /*ひび割れのない単鉄筋コンクリート */
    int mm=3;
    double sm[Mmax+1][Mmax+2];
    double Ac,Bc,C1,C2,Cg,Psa;
    int i,j;

    STIF(STRES,Ecc,ncc,&Ac,&Bc);

    for(i=0;i<=mm;i++){
        for(j=0;j<=mm+1;j++){
            sm[i][j]=0.0;
    }}
    sm[0][0]=ar;
    sm[0][1]=1.0/ar;
    sm[0][2]=-ar*ar/Ess/taa;
    sm[1][0]=Ac+Bc;
    sm[1][1]=-(Ac-Bc)/ar/ar;
    sm[1][2]=-1.0;
    sm[2][0]=br;
    sm[2][1]=1.0/br;
    sm[2][3]=-1.0/br;
    sm[3][0]=Ac+Bc;
    sm[3][1]=-(Ac-Bc)/br/br;
    sm[3][3]=Egg/(1.0+ngg)/br/br;
    sm[1][4]=-PP;
    MATGJ(mm,sm);
    C1=sm[0][4];
    C2=sm[1][4];
    Psa=sm[2][4];
    Cg=sm[3][4];
    *ssa=Psa*ar/taa;
    *ssb=0.0;
    *src=(Ac+Bc)*C1-(Ac-Bc)*C2/ar/ar;
    *stc=(Ac+Bc)*C1+(Ac-Bc)*C2/ar/ar;
    *uaa=Psa*ar*ar/Ess/taa;
    *ubb=C1*br+C2/br;

}
/*----------------------------------------------------------------------------*/

void CALC2(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb)
{
    /*ひび割れのある単鉄筋コンクリート */
    int mm=3;
    double sm[Mmax+1][Mmax+2];
    double C3,C4,Cg,Psa;
    int i,j;

    for(i=0;i<=mm;i++){
        for(j=0;j<=mm+1;j++){
            sm[i][j]=0.0;
    }}
    sm[0][0]=log(ar)/Ecc;
    sm[0][1]=1.0;
    sm[0][2]=-ar*ar/Ess/taa;
    sm[1][0]=1.0/ar;
    sm[1][2]=-1.0;
    sm[2][0]=log(br)/Ecc;
    sm[2][1]=1.0;
    sm[2][3]=-1.0/br;
    sm[3][0]=1.0/br;
    sm[3][3]=Egg/(1.0+ngg)/br/br;
    sm[1][4]=-PP;
    MATGJ(mm,sm);
    C3=sm[0][4];
    C4=sm[1][4];
    Psa=sm[2][4];
    Cg=sm[3][4];
    *ssa=Psa*ar/taa;
    *ssb=0.0;
    *src=C3/ar;
    *stc=0.0;
    *uaa=Psa*ar*ar/Ess/taa;
    *ubb=C3/Ecc*log(br)+C4;

}
/*----------------------------------------------------------------------------*/

void CALC3(int STRES,double PP,double ar,double br,
    double Ecc,double ncc,double Ess,double taa,double tbb,double Egg,double ngg,
    double *ssa,double *ssb,double *src,double *stc,double *uaa,double *ubb)
{
    /*ひび割れのある複鉄筋コンクリート */
    int mm=4;
    double sm[Mmax+1][Mmax+2];
    double C3,C4,Cg,Psa,Psb;
    int i,j;
    
    if(0.0<tbb){
        for(i=0;i<=mm;i++){
            for(j=0;j<=mm+1;j++){
                sm[i][j]=0.0;
        }}
        sm[0][0]=log(ar)/Ecc;
        sm[0][1]=1.0;
        sm[0][2]=-ar*ar/Ess/taa;
        sm[1][0]=1.0/ar;
        sm[1][2]=-1.0;
        sm[2][0]=log(br)/Ecc;
        sm[2][1]=1.0;
        sm[2][3]=-br*br/Ess/tbb;
        sm[3][3]=br*br/Ess/tbb;
        sm[3][4]=-1.0/br;
        sm[4][0]=1.0/br;
        sm[4][3]=1.0;
        sm[4][4]=Egg/(1.0+ngg)/br/br;
        sm[1][5]=-PP;
        MATGJ(mm,sm);
        C3=sm[0][5];
        C4=sm[1][5];
        Psa=sm[2][5];
        Psb=sm[3][5];
        Cg=sm[4][5];
        *ssa=Psa*ar/taa;
        *ssb=Psb*br/tbb;
        *src=C3/ar;
        *stc=0.0;
        *uaa=Psa*ar*ar/Ess/taa;
        *ubb=C3/Ecc*log(br)+C4;
    }else{
        *ssa=0.0;
        *ssb=0.0;
        *src=0.0;
        *stc=0.0;
        *uaa=0.0;
        *ubb=0.0;
    }

}
/*----------------------------------------------------------------------------*/

void MATGJ(int n,double a[Mmax+1][Mmax+2])
{
    /* Gauss-Jordan法による連立一次方程式の解法 */
    int i,j,k,s;
    double p,d,max,dumy;

    for(k=0;k<=n;k++){
        max=0.0;
        s=k;
        for(j=k;j<=n;j++){
            if(fabs(a[j][k])>max){
                max=fabs(a[j][k]);
                s=j;
            }
        }
        for(j=0;j<=n+1;j++){
            dumy=a[k][j];
            a[k][j]=a[s][j];
            a[s][j]=dumy;
        }
        p=a[k][k];
        for(j=k;j<=n+1;j++){
            a[k][j]=a[k][j]/p;
        }
        for(i=0;i<=n;i++){
            if(i!=k){
                d=a[i][k];
                for(j=k;j<=n+1;j++){
                    a[i][j]=a[i][j]-d*a[k][j];
                }
            }
        }
    }
}
/*-------------------------------------------------------------------*/