/*===========================*/
/* gccCPTT.c                 */
/*===========================*/
/* Stress analysis for Circular Pressure Tunnel */
/*     using multiple cylindrical theorem */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 12
#define NDmax 100

char strcom[256];       /* Comment */
int    nd;              /* Numbers of calculation */
int    IE_inp[NDmax+1]; /* 0:Inner pressure, 1:External pressure */
double PP_inp[NDmax+1]; /* Water pressure */
double TT_inp[NDmax+1]; /* Temperature change (+:increse of temperature) */
double aa_inp[NDmax+1]; /* Internal radius of pressure tunnel */
double bb_inp[NDmax+1]; /* Excavation radius */
double rr_inp[NDmax+1]; /* Model edge radius */
double cc_inp[NDmax+1]; /* Cover for reinforcement */
double ta_inp[NDmax+1]; /* Equivalent thickness of inner reinforcement */
double tb_inp[NDmax+1]; /* Equivalent thickness of outer reinforcement */
double Ec_inp[NDmax+1]; /* Elastic modulus of concrete */
double nc_inp[NDmax+1]; /* Poisson ratio of concrete */
double ac_inp[NDmax+1]; /* Thermal expansion coefficient of concrete */
double Es_inp[NDmax+1]; /* Elastic modulus of reinforcement */
double ns_inp[NDmax+1]; /* Poisson ratio of reinforcement */
double as_inp[NDmax+1]; /* Thermal expansion coefficient of reinforcement */
double Eg_inp[NDmax+1]; /* Elastic modulus of bed rock */
double ng_inp[NDmax+1]; /* Poisson ratio of bed rock */

void main_part(char fnameR[],char fnameW[]);
void main_sng_pin(char fnameW[],int k);
void main_dbl_pin(char fnameW[],int k);
void main_sng_pex(char fnameW[],int k);
void main_dbl_pex(char fnameW[],int k);
void ELAS(double C1,double C2,double xx,double xa,double E,double nu,double alp,double TT,
          double *u,double *sigr,double *sigt);
void CON(double C1,double C2,double xx,double xa,double E,double alp,double TT,double *u,double *sigr);
void MATGJ(int n,double a[Nmax][Nmax+1]);
/*------------------------------------------------------------------------------*/
int main(int argc,char *argv[])
{
    char fnameR[50]; /* Input file name */
    char fnameW[50]; /* Output file name */

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

    main_part(fnameR,fnameW);

    return 0;
}
/*------------------------------------------------------------------------------*/
void main_part(char fnameR[],char fnameW[])
{
    FILE *fin,*fout;
    char dat[1000];
    char sv[20];
    int i,k;

    /* Data input */
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);strcpy(strcom,strtok(dat,","));
    k=0;
    while(fgets(dat,sizeof dat,fin)!=NULL){
        if((strstr(dat,"#"))==NULL){
            IE_inp[k]=atoi(strtok(dat,","));
            PP_inp[k]=atof(strtok(NULL,","));
            TT_inp[k]=atof(strtok(NULL,","));
            aa_inp[k]=atof(strtok(NULL,","));
            bb_inp[k]=atof(strtok(NULL,","));
            rr_inp[k]=atof(strtok(NULL,","));
            cc_inp[k]=atof(strtok(NULL,","));
            ta_inp[k]=atof(strtok(NULL,","));
            tb_inp[k]=atof(strtok(NULL,","));
            Ec_inp[k]=atof(strtok(NULL,","));
            nc_inp[k]=atof(strtok(NULL,","));
            ac_inp[k]=atof(strtok(NULL,","));
            Es_inp[k]=atof(strtok(NULL,","));
            ns_inp[k]=atof(strtok(NULL,","));
            as_inp[k]=atof(strtok(NULL,","));
            Eg_inp[k]=atof(strtok(NULL,","));
            ng_inp[k]=atof(strtok(NULL,","));
            k=k+1;
        }
    }
    fclose(fin);
    nd=k-1;

    /* Making output file */
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s",strcom);
    fprintf(fout,"*Input data\n");
    fprintf(fout,"k,IE,PP,TT,aa,bb,rr,cc,ta,tb,Ec,nc,ac,Es,ns,as,Eg,ng\n");
    for(i=0;i<=nd;i++){
        sprintf(sv,"%d",i);strcpy(dat,sv);
        sprintf(sv,"%d",IE_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",PP_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",TT_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",aa_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",bb_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",rr_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",cc_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",ta_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",tb_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",Ec_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",nc_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",ac_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",Es_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",ns_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",as_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",Eg_inp[i]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%.3f",ng_inp[i]);strcat(dat,",");strcat(dat,sv);
        strcat(dat,"\n");
        fputs(dat,fout);
    }
    fprintf(fout,"*Output data\n");
    /* k     : Sequence                                                     */
    /* IE    : Inner pressure (=0) or Outer pressure (=1)                   */
    /* Eg    : Water pressure                                               */
    /* dT    : Temperature change                                           */
    /* sr_c  : Stress in the r-direction at inner surface concrete             */
    /* st_c  : Stress in the t-direction at inner surface concrete             */
    /* sr_si : Stress in the r-direction of inner reinforcement                */
    /* st_si : Stress in the t-direction of inner reinforcement                */
    /* sr_so : Stress in the r-direction of outer reinforcement                */
    /* st_so : Stress in the t-direction of outer reinforcement                */
    /* sr_g  : Stress in the r-direction at outer surface concrete or bed rock */
    /* st_g  : Stress in the t-direction at outer surface concrete or bed rock */
    /* ua    : Radial displacement of inner surface concrete                */
    /* ub    : Radial displacement of outer surface concrete                */
    fprintf(fout,"k,IE,Eg,dT,sr_c,st_c,sr_si,st_si,sr_so,st_so,sr_g,st_g,ua,ub\n");
    fclose(fout);

    /* Execution of calculation */
    for(k=0;k<=nd;k++){
        switch(IE_inp[k]){
            case 0:
                if(tb_inp[k]<=0)main_sng_pin(fnameW,k);
                if(0<tb_inp[k]) main_dbl_pin(fnameW,k);
            break;
            case 1:
                if(tb_inp[k]<=0)main_sng_pex(fnameW,k);
                if(0<tb_inp[k]) main_dbl_pex(fnameW,k);
            break;
        }
    }

}
/*------------------------------------------------------------------------------*/
void main_sng_pin(char fnameW[],int k)
{
    FILE *fout;
    int n=7;
    int i,j;
    double a[Nmax][Nmax+1];
    double ra,r1,r2,rb,rr;
    double Eg,ng,Es,ns,alps,Ec,alpc;
    double P0,DT;
    double C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt;
    double ug,sigrg,sigtg;
    double us2,sigrs2,sigts2;
    double us1,sigrs1,sigts1;
    double uc,sigrc;

    ra  =aa_inp[k];
    r1  =aa_inp[k]+cc_inp[k];
    r2  =r1+ta_inp[k];
    rb  =bb_inp[k];
    rr  =rr_inp[k];
    Eg  =Eg_inp[k];
    ng  =ng_inp[k];
    Es  =Es_inp[k];
    ns  =ns_inp[k];
    alps=as_inp[k];
    Ec  =Ec_inp[k];
    alpc=ac_inp[k];
    P0  =PP_inp[k]; /* Positive : internal pressure */
    DT  =TT_inp[k];

    for(i=0;i<=n;i++){
        for(j=0;j<=n+1;j++){
            a[i][j]=0.0;
    }}

    a[0][0]=rr                      ;a[0][1]= 1.0/rr;
    a[1][0]=Eg/(1.0+ng)/(1.0-2.0*ng);a[1][1]=-Eg/(1.0+ng)/rb/rb;a[1][2]= 0.0                     ;a[1][3]=-Ec/rb;
    a[2][0]=rb                      ;a[2][1]= 1.0/rb           ;a[2][2]=-1.0                     ;a[2][3]=-log(rb);
    a[3][2]=0.0                     ;a[3][3]= Ec/r2            ;a[3][4]=-Es/(1.0+ns)/(1.0-2.0*ns);a[3][5]= Es/(1.0+ns)/r2/r2;
    a[4][2]=1.0                     ;a[4][3]= log(r2)          ;a[4][4]=-r2                      ;a[4][5]=-1.0/r2;
    a[5][4]=Es/(1.0+ns)/(1.0-2.0*ns);a[5][5]=-Es/(1.0+ns)/r1/r1;a[5][6]= 0.0                     ;a[5][7]=-Ec/r1;
    a[6][4]=r1                      ;a[6][5]= 1.0/r1           ;a[6][6]=-1.0                     ;a[6][7]=-log(r1);
    a[7][6]=0.0                     ;a[7][7]= Ec/ra;

    a[0][8]= 0.0;
    a[1][8]= 0.0;
    a[2][8]= alpc*DT*(rb-r2);
    a[3][8]=-Es*alps*DT/(1.0-ns)*(r2*r2-r1*r1)/2.0/r2/r2;
    a[4][8]= (1.0+ns)/(1.0-ns)*alps*DT*(r2*r2-r1*r1)/2.0/r2;
    a[5][8]= 0.0;
    a[6][8]= alpc*DT*(r1-ra);
    a[7][8]=-P0;

    MATGJ(n,a);

    C1=a[0][8];C2=a[1][8];xx=rb;xa=rb;E=Eg;nu=ng;alp=0.0;TT=0.0;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    ug=u;sigrg=sigr;sigtg=sigt;

    C1=a[4][8];C2=a[5][8];xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us2=u;sigrs2=sigr;sigts2=sigt;

    C1=a[4][8];C2=a[5][8];xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us1=u;sigrs1=sigr;sigts1=sigt;

    C1=a[6][8];C2=a[7][8];xx=ra;xa=ra;E=Ec;alp=alpc;TT=DT;
    CON(C1,C2,xx,xa,E,alp,DT,&u,&sigr);
    uc=u;sigrc=sigr;

    fout=fopen(fnameW,"a");
    fprintf(fout,"%d,%d,%.0f,%.0f,%.3f,%.3f,%.3f,%.3f,---,---,%.3f,%.3f,%.3f,%.3f\n",
        k,IE_inp[k],Eg,DT,sigrc,0.0,0.5*(sigrs1+sigrs2),0.5*(sigts1+sigts2),sigrg,sigtg,uc,ug);
    fclose(fout);

    printf("===Single: Internal pressure===\n");
    printf("%.3f %.3f ---\n",uc,sigrc);
    printf("%.3f %.3f %.3f\n",us1,sigrs1,sigts1);
    printf("%.3f %.3f %.3f\n",us2,sigrs2,sigts2);
    printf("%.3f %.3f %.3f\n",ug,sigrg,sigtg);

}
/*------------------------------------------------------------------------------*/
void main_dbl_pin(char fnameW[],int k)
{
    FILE *fout;
    int n=11;
    int i,j;
    double a[Nmax][Nmax+1];
    double ra,r1,r2,r3,r4,rb,rr;
    double Eg,ng,Es,ns,alps,Ec,alpc;
    double P0,DT;
    double C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt;
    double ug,sigrg,sigtg;
    double us4,sigrs4,sigts4;
    double us3,sigrs3,sigts3;
    double us2,sigrs2,sigts2;
    double us1,sigrs1,sigts1;
    double uc,sigrc;

    int nnn;

    ra  =aa_inp[k];
    r1  =ra+cc_inp[k];
    r2  =r1+ta_inp[k];
    rb  =bb_inp[k];
    r4  =rb-cc_inp[k];
    r3  =r4-tb_inp[k];
    rr  =rr_inp[k];
    Eg  =Eg_inp[k];
    ng  =ng_inp[k];
    Es  =Es_inp[k];
    ns  =ns_inp[k];
    alps=as_inp[k];
    Ec  =Ec_inp[k];
    alpc=ac_inp[k];
    P0  =PP_inp[k]; /* Positive : internal pressure */
    DT  =TT_inp[k];

    for(i=0;i<=n;i++){
        for(j=0;j<=n+1;j++){
            a[i][j]=0.0;
    }}

    a[0][0]  =rr                      ;a[0][1]  = 1.0/rr;
    a[1][0]  =Eg/(1.0+ng)/(1.0-2.0*ng);a[1][1]  =-Eg/(1.0+ng)/rb/rb ;a[1][2]  = 0.0                     ;a[1][3]  =-Ec/rb;
    a[2][0]  =rb                      ;a[2][1]  = 1.0/rb            ;a[2][2]  =-1.0                     ;a[2][3]  =-log(rb);
    a[3][2]  =0.0                     ;a[3][3]  = Ec/r4             ;a[3][4]  =-Es/(1.0+ns)/(1.0-2.0*ns);a[3][5]  = Es/(1.0+ns)/r4/r4;
    a[4][2]  =1.0                     ;a[4][3]  = log(r4)           ;a[4][4]  =-r4                      ;a[4][5]  =-1.0/r4;
    a[5][4]  =Es/(1.0+ns)/(1.0-2.0*ns);a[5][5]  =-Es/(1.0+ns)/r3/r3 ;a[5][6]  = 0.0                     ;a[5][7]  =-Ec/r3;
    a[6][4]  =r3                      ;a[6][5]  = 1.0/r3            ;a[6][6]  =-1.0                     ;a[6][7]  =-log(r3);
    a[7][6]  =0.0                     ;a[7][7]  = Ec/r2             ;a[7][8]  =-Es/(1.0+ns)/(1.0-2.0*ns);a[7][9]  = Es/(1.0+ns)/r2/r2;
    a[8][6]  =1.0                     ;a[8][7]  = log(r2)           ;a[8][8]  =-r2                      ;a[8][9]  =-1.0/r2;
    a[9][8]  =Es/(1.0+ns)/(1.0-2.0*ns);a[9][9]  =-Es/(1.0+ns)/r1/r1 ;a[9][10] = 0.0                     ;a[9][11] =-Ec/r1;
    a[10][8] =r1                      ;a[10][9] = 1.0/r1            ;a[10][10]=-1.0                     ;a[10][11]=-log(r1);
    a[11][10]=0.0                     ;a[11][11]= Ec/ra;

    a[0][12] = 0.0;
    a[1][12] = 0.0;
    a[2][12] = alpc*DT*(rb-r2);
    a[3][12] =-Es*alps*DT/(1.0-ns)*(r4*r4-r3*r3)/2.0/r4/r4;
    a[4][12] = (1.0+ns)/(1.0-ns)*alps*DT*(r4*r4-r3*r3)/2.0/r4;
    a[5][12] = 0.0;
    a[6][12] = alpc*DT*(r3-r2);
    a[7][12] =-Es*alps*DT/(1.0-ns)*(r2*r2-r1*r1)/2.0/r2/r2;
    a[8][12] = (1.0+ns)/(1.0-ns)*alps*DT*(r2*r2-r1*r1)/2.0/r2;
    a[9][12] = 0.0;
    a[10][12]= alpc*DT*(r1-ra);
    a[11][12]=-P0;

    MATGJ(n,a);

    C1=a[0][12];C2=a[1][12];xx=rb;xa=rb;E=Eg;nu=ng;alp=0.0;TT=0.0;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    ug=u;sigrg=sigr;sigtg=sigt;

    C1=a[4][12];C2=a[5][12];xx=r4;xa=r3;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us4=u;sigrs4=sigr;sigts4=sigt;

    C1=a[4][12];C2=a[5][12];xx=r3;xa=r3;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us3=u;sigrs3=sigr;sigts3=sigt;

    C1=a[8][12];C2=a[9][12];xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us2=u;sigrs2=sigr;sigts2=sigt;

    C1=a[8][12];C2=a[9][12];xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us1=u;sigrs1=sigr;sigts1=sigt;

    C1=a[10][12];C2=a[11][12];xx=ra;xa=ra;E=Ec;alp=alpc;TT=DT;
    CON(C1,C2,xx,xa,E,alp,DT,&u,&sigr);
    uc=u;sigrc=sigr;

    fout=fopen(fnameW,"a");
    fprintf(fout,"%d,%d,%.0f,%.0f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n",
        k,IE_inp[k],Eg,DT,sigrc,0.0,0.5*(sigrs1+sigrs2),0.5*(sigts1+sigts2),0.5*(sigrs3+sigrs4),0.5*(sigts3+sigts4),sigrg,sigtg,uc,ug);
    fclose(fout);

    printf("===Double: Internal pressure===\n");
    printf("%.3f %.3f ---\n",uc,sigrc);
    printf("%.3f %.3f %.3f\n",us1,sigrs1,sigts1);
    printf("%.3f %.3f %.3f\n",us2,sigrs2,sigts2);
    printf("%.3f %.3f %.3f\n",us3,sigrs3,sigts3);
    printf("%.3f %.3f %.3f\n",us4,sigrs4,sigts4);
    printf("%.3f %.3f %.3f\n",ug,sigrg,sigtg);

}
/*------------------------------------------------------------------------------*/
void main_sng_pex(char fnameW[],int k)
{
    FILE *fout;
    int n=5;
    int i,j;
    double a[Nmax][Nmax+1];
    double ra,r1,r2,rb,rr;
    double Eg=0.0,Es,ns,alps,Ec,nc,alpc;
    double Pe,DT;
    double C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt;
    double ug,sigrg,sigtg;
    double us2,sigrs2,sigts2;
    double us1,sigrs1,sigts1;
    double uc,sigrc,sigtc;

    ra  =aa_inp[k];
    r1  =aa_inp[k]+cc_inp[k];
    r2  =r1+ta_inp[k];
    rb  =bb_inp[k];
    Es  =Es_inp[k];
    ns  =ns_inp[k];
    alps=as_inp[k];
    Ec  =Ec_inp[k];
    nc  =nc_inp[k];
    alpc=ac_inp[k];
    Pe  =PP_inp[k]; /* Positive : outer pressure */
    DT  =TT_inp[k];

    for(i=0;i<=n;i++){
        for(j=0;j<=n+1;j++){
            a[i][j]=0.0;
    }}

    a[0][0]=Ec/(1.0+nc)/(1.0-2.0*nc);a[0][1]=-Ec/(1.0+nc)/rb/rb;
    a[1][0]=Ec/(1.0+nc)/(1.0-2.0*nc);a[1][1]=-Ec/(1.0+nc)/r2/r2;a[1][2]=-Es/(1.0+ns)/(1.0-2.0*ns);a[1][3]=Es/(1.0+ns)/r2/r2;
    a[2][0]=r2                      ;a[2][1]= 1.0/r2           ;a[2][2]=-r2                      ;a[2][3]=-1.0/r2;
    a[3][2]=Es/(1.0+ns)/(1.0-2.0*ns);a[3][3]=-Es/(1.0+ns)/r1/r1;a[3][4]=-Ec/(1.0+nc)/(1.0-2.0*nc);a[3][5]=Ec/(1.0+nc)/r1/r1;
    a[4][2]=r1                      ;a[4][3]= 1.0/r1           ;a[4][4]=-r1                      ;a[4][5]=-1.0/r1;
    a[5][4]=Ec/(1.0+nc)/(1.0-2.0*nc);a[5][5]=-Ec/(1.0+nc)/ra/ra;

    a[0][6]=-Pe;
    a[1][6]=-Es*alps*DT/(1.0-ns)*(r2*r2-r1*r1)/2.0/r2/r2;
    a[2][6]= (1.0+ns)/(1.0-ns)*alps*DT*(r2*r2-r1*r1)/2.0/r2;
    a[3][6]=-Ec*alpc*DT/(1.0-nc)*(r1*r1-ra*ra)/2.0/r1/r1;
    a[4][6]= (1.0+nc)/(1.0-nc)*alpc*DT*(r1*r1-ra*ra)/2.0/r1;
    a[5][6]= 0.0;

    MATGJ(n,a);

    C1=a[0][6];C2=a[1][6];xx=rb;xa=r2;E=Ec;nu=nc;alp=alpc;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    ug=u;sigrg=sigr;sigtg=sigt;

    C1=a[2][6];C2=a[3][6];xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us2=u;sigrs2=sigr;sigts2=sigt;

    C1=a[2][6];C2=a[3][6];xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us1=u;sigrs1=sigr;sigts1=sigt;

    C1=a[4][6];C2=a[5][6];xx=ra;xa=ra;E=Ec;nu=nc;alp=alpc;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    uc=u;sigrc=sigr;sigtc=sigt;

    fout=fopen(fnameW,"a");
    fprintf(fout,"%d,%d,%.0f,%.0f,%.3f,%.3f,%.3f,%.3f,---,---,%.3f,%.3f,%.3f,%.3f\n",
        k,IE_inp[k],Eg,DT,sigrc,sigtc,0.5*(sigrs1+sigrs2),0.5*(sigts1+sigts2),sigrg,sigtg,uc,ug);
    fclose(fout);

    printf("===Single: External pressure===\n");
    printf("%.3f %.3f %.3f\n",uc,sigrc,sigtc);
    printf("%.3f %.3f %.3f\n",us1,sigrs1,sigts1);
    printf("%.3f %.3f %.3f\n",us2,sigrs2,sigts2);
    printf("%.3f %.3f %.3f\n",ug,sigrg,sigtg);

}
/*------------------------------------------------------------------------------*/
void main_dbl_pex(char fnameW[],int k)
{
    FILE *fout;
    int n=9;
    int i,j;
    double a[Nmax][Nmax+1];
    double ra,r1,r2,r3,r4,rb,rr;
    double Eg=0.0,Es,ns,alps,Ec,nc,alpc;
    double Pe,DT;
    double C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt;
    double ug,sigrg,sigtg;
    double us4,sigrs4,sigts4;
    double us3,sigrs3,sigts3;
    double us2,sigrs2,sigts2;
    double us1,sigrs1,sigts1;
    double uc,sigrc,sigtc;

    int nnn;

    ra  =aa_inp[k];
    r1  =ra+cc_inp[k];
    r2  =r1+ta_inp[k];
    rb  =bb_inp[k];
    r4  =rb-cc_inp[k];
    r3  =r4-tb_inp[k];
    Es  =Es_inp[k];
    ns  =ns_inp[k];
    alps=as_inp[k];
    Ec  =Ec_inp[k];
    nc  =nc_inp[k];
    alpc=ac_inp[k];
    Pe  =PP_inp[k]; /* Positive : outer pressure */
    DT  =TT_inp[k];

    for(i=0;i<=n;i++){
        for(j=0;j<=n+1;j++){
            a[i][j]=0.0;
    }}

    a[0][0]=Ec/(1.0+nc)/(1.0-2.0*nc);a[0][1]=-Ec/(1.0+nc)/rb/rb;
    a[1][0]=Ec/(1.0+nc)/(1.0-2.0*nc);a[1][1]=-Ec/(1.0+nc)/r4/r4;a[1][2]=-Es/(1.0+ns)/(1.0-2.0*ns);a[1][3]=Es/(1.0+ns)/r4/r4;
    a[2][0]=r4                      ;a[2][1]= 1.0/r4           ;a[2][2]=-r4                      ;a[2][3]=-1.0/r4;
    a[3][2]=Es/(1.0+ns)/(1.0-2.0*ns);a[3][3]=-Es/(1.0+ns)/r3/r3;a[3][4]=-Ec/(1.0+nc)/(1.0-2.0*nc);a[3][5]=Ec/(1.0+nc)/r3/r3;
    a[4][2]=r3                      ;a[4][3]= 1.0/r3           ;a[4][4]=-r3                      ;a[4][5]=-1.0/r3;
    a[5][4]=Ec/(1.0+nc)/(1.0-2.0*nc);a[5][5]=-Ec/(1.0+nc)/r2/r2;a[5][6]=-Es/(1.0+ns)/(1.0-2.0*ns);a[5][7]=Es/(1.0+ns)/r2/r2;
    a[6][4]=r2                      ;a[6][5]= 1.0/r2           ;a[6][6]=-r2                      ;a[6][7]=-1.0/r2;
    a[7][6]=Es/(1.0+ns)/(1.0-2.0*ns);a[7][7]=-Es/(1.0+ns)/r1/r1;a[7][8]=-Ec/(1.0+nc)/(1.0-2.0*nc);a[7][9]=Ec/(1.0+nc)/r1/r1;
    a[8][6]=r1                      ;a[8][7]= 1.0/r1           ;a[8][8]=-r1                      ;a[8][9]=-1.0/r1;
    a[9][8]=Ec/(1.0+nc)/(1.0-2.0*nc);a[9][9]=-Ec/(1.0+nc)/ra/ra;

    a[0][10]=-Pe;
    a[1][10]=-Es*alps*DT/(1.0-ns)*(r4*r4-r3*r3)/2.0/r4/r4;
    a[2][10]= (1.0+ns)/(1.0-ns)*alps*DT*(r4*r4-r3*r3)/2.0/r4;
    a[3][10]=-Ec*alpc*DT/(1.0-nc)*(r3*r3-r2*r2)/2.0/r3/r3;
    a[4][10]= (1.0+nc)/(1.0-nc)*alpc*DT*(r3*r3-r2*r2)/2.0/r3;
    a[5][10]=-Es*alps*DT/(1.0-ns)*(r2*r2-r1*r1)/2.0/r2/r2;
    a[6][10]= (1.0+ns)/(1.0-ns)*alps*DT*(r2*r2-r1*r1)/2.0/r2;
    a[7][10]=-Ec*alpc*DT/(1.0-nc)*(r1*r1-ra*ra)/2.0/r1/r1;
    a[8][10]= (1.0+nc)/(1.0-nc)*alpc*DT*(r1*r1-ra*ra)/2.0/r1;
    a[9][10]= 0.0;

    MATGJ(n,a);

    C1=a[0][10];C2=a[1][10];xx=rb;xa=r4;E=Ec;nu=nc;alp=alpc;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    ug=u;sigrg=sigr;sigtg=sigt;

    C1=a[2][10];C2=a[3][10];xx=r4;xa=r3;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us4=u;sigrs4=sigr;sigts4=sigt;

    C1=a[2][10];C2=a[3][10];xx=r3;xa=r3;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us3=u;sigrs3=sigr;sigts3=sigt;

    C1=a[6][10];C2=a[7][10];xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us2=u;sigrs2=sigr;sigts2=sigt;

    C1=a[6][10];C2=a[7][10];xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    us1=u;sigrs1=sigr;sigts1=sigt;

    C1=a[8][10];C2=a[9][10];xx=ra;xa=ra;E=Ec;nu=nc;alp=alpc;TT=DT;
    ELAS(C1,C2,xx,xa,E,nu,alp,TT,&u,&sigr,&sigt);
    uc=u;sigrc=sigr;sigtc=sigt;

    fout=fopen(fnameW,"a");
    fprintf(fout,"%d,%d,%.0f,%.0f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n",
        k,IE_inp[k],Eg,DT,sigrc,sigtc,0.5*(sigrs1+sigrs2),0.5*(sigts1+sigts2),0.5*(sigrs3+sigrs4),0.5*(sigts3+sigts4),sigrg,sigtg,uc,ug);
    fclose(fout);

    printf("===Double: External pressure===\n");
    printf("%.3f %.3f %.3f\n",uc,sigrc,sigtc);
    printf("%.3f %.3f %.3f\n",us1,sigrs1,sigts1);
    printf("%.3f %.3f %.3f\n",us2,sigrs2,sigts2);
    printf("%.3f %.3f %.3f\n",us3,sigrs3,sigts3);
    printf("%.3f %.3f %.3f\n",us4,sigrs4,sigts4);
    printf("%.3f %.3f %.3f\n",ug,sigrg,sigtg);

}
/*------------------------------------------------------------------------------*/
void ELAS(double C1,double C2,double xx,double xa,double E,double nu,double alp,double TT,
          double *u,double *sigr,double *sigt)
{
    *u=(1.0+nu)/(1.0-nu)*alp*TT*(xx*xx-xa*xa)/2.0/xx+C1*xx+C2/xx;
    *sigr=-E*alp*TT/(1.0-nu)*(xx*xx-xa*xa)/2.0/xx/xx+E/(1.0+nu)/(1.0-2.0*nu)*C1-E/(1.0+nu)*C2/xx/xx;
    *sigt=-E*alp*TT/(1.0-nu)*(xx*xx+xa*xa)/2.0/xx/xx+E/(1.0+nu)/(1.0-2.0*nu)*C1+E/(1.0+nu)*C2/xx/xx;
}
/*------------------------------------------------------------------------------*/
void CON(double C1,double C2,double xx,double xa,double E,double alp,double TT,double *u,double *sigr)
{
    *u=C1+C2*log(xx)+alp*TT*(xx-xa);
    *sigr=E*C2/xx;
}
/*------------------------------------------------------------------------------*/
void MATGJ(int n,double a[Nmax][Nmax+1])
{
    /* 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];
                }
            }
        }
    }
}
/*------------------------------------------------------------------------------*/