/************************************************/
/* gccFIG_PSCPT.c (埋設鉄管・RC圧力水路) */
/************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define pi M_PI
#define Nmax 200

void main_part(int iii,char fnameW[50],char fnameF[50],double rin,double rout,double t0,
               int n,double R,double L,double a,
               double r1,double r2,double t1,double t2,double dr1,double dr2);
void STR_PS(FILE *fp,int iii,double rin,double rout,double t0);
void STR_RC(FILE *fp,int iii,double rin,double rout,double t0);
void AGR(FILE *fp,double r1,double r2,double t1,double t2,double dr1,double dr2);
void SYMB_CRACK(FILE *fp,int n,double R,double L,double a);
void SYMB_ROCK(FILE *fp,double rout);
void CAL_TR(double phi,double R,double x1,double y1,double *co_r,double *co_t);
void CAL_XY(double r1,double r2,double t1,double t2,double *xx,double *yy);
void CAL_XYARRAY(double r1,double r2,double t1,double t2,double dr1,double dr2,double x[Nmax],double y[Nmax],int *MG,int *MT);
/*---------------------------------------------------------------------------*/
int main()
{
    char fnameW[50];
    char fnameF[50];
    int iii,n;
    double rin,rout,t0;
    double phi,R,L,a;
    double r1,r2,t1,t2,dr1,dr2;

    /* ひび割れ */
    n=4;
    a=0.1;

    /* 骨材 */
    dr1=1.0;
    dr2=0.5;

    for(iii=0;iii<=4;iii++){
        if(iii==0){
            /* 構造形状(埋設鉄管) */
            rin=4.0;
            rout=7.0;
            t0=0.5;
            /* 骨材描画範囲 */
            t1=45;
            t2=225;
        }else{
            /* 構造形状(RC水路) */
            rin=4.5;
            rout=7.0;
            t0=0.2;
            /* 骨材描画範囲 */
            t1=-45;
            t2=225;
        }
        switch(iii){
            case 0: /* Penstock */
                strcpy(fnameW,"gpl_eps_model_PS.txt");
                strcpy(fnameF,"fig_eps_model_PS.eps");
                R=rin+t0;      /* ひび割れ */
                L=rout-t0-rin; /* ひび割れ */
                r1=rin+t0+0.2; /* 骨材 */
                r2=rout-0.2;   /* 骨材 */
                break;
            case 1: /* No-reinforcement no-crack */
                strcpy(fnameW,"gpl_eps_model_CPT1.txt");
                strcpy(fnameF,"fig_eps_model_CPT1.eps");
                r1=rin+0.2;  /* 骨材 */
                r2=rout-0.2; /* 骨材 */
                break;
            case 2: /* Single reinforcement no-crack */
                strcpy(fnameW,"gpl_eps_model_CPT2.txt");
                strcpy(fnameF,"fig_eps_model_CPT2.eps");
                r1=rin+t0+0.2; /* 骨材 */
                r2=rout-0.2;   /* 骨材 */
                break;
            case 3: /* Single reinforcements with cracks */
                strcpy(fnameW,"gpl_eps_model_CPT3.txt");
                strcpy(fnameF,"fig_eps_model_CPT3.eps");
                R=rin+t0;       /* ひび割れ */
                L=rout-t0-rin;  /* ひび割れ */
                r1=rin+t0+0.2;  /* 骨材 */
                r2=rout-0.2;    /* 骨材 */
                break;
            case 4: /* Double reinforcements with cracks */
                strcpy(fnameW,"gpl_eps_model_CPT4.txt");
                strcpy(fnameF,"fig_eps_model_CPT4.eps");
                R=rin+t0;           /* ひび割れ */
                L=rout-t0-t0-rin;  /* ひび割れ */
                r1=rin+t0+0.2;     /* 骨材 */
                r2=rout-t0-0.2;    /* 骨材 */
                break;
        }
        main_part(iii,fnameW,fnameF,rin,rout,t0,n,R,L,a,r1,r2,t1,t2,dr1,dr2);
    }

}
/*---------------------------------------------------------------------------*/
void main_part(int iii,char fnameW[50],char fnameF[50],double rin,double rout,double t0,
               int n,double R,double L,double a,
               double r1,double r2,double t1,double t2,double dr1,double dr2)
{
    FILE *fp;

    fp=fopen(fnameW,"w");
    if(iii==0){
        fprintf(fp,"set terminal postscript eps enhanced font \"Helvetica\" 20\n");
    }else{
        fprintf(fp,"set terminal postscript eps enhanced font \"Helvetica\" 26\n");
    }
    fprintf(fp,"set output \"%s\"\n",fnameF);
    fprintf(fp,"#\n");
    fprintf(fp,"set multiplot\n");
    fprintf(fp,"set origin 0,0\n");
    fprintf(fp,"#\n");
    fprintf(fp,"set polar\n");
    fprintf(fp,"set angles degree\n");
    fprintf(fp,"set size ratio -1\n");
    fprintf(fp,"set samples 360\n");
    fprintf(fp,"#\n");
    fprintf(fp,"set xrange [-10:10]\n");
    fprintf(fp,"set yrange [-10:10]\n");
    fprintf(fp,"set trange [0:360]\n");
    fprintf(fp,"set border 0\n");
    fprintf(fp,"set xtics 0,0,0\n");
    fprintf(fp,"set ytics 0,0,0\n");
    fprintf(fp,"set xtics nomirror\n");
    fprintf(fp,"set ytics nomirror\n");
    fprintf(fp,"#\n");
    if(iii==0){
        STR_PS(fp,iii,rin,rout,t0); /*埋設鉄管-文字描画*/
    }else{
        STR_RC(fp,iii,rin,rout,t0); /*RC水路-文字描画*/
    }
    /*骨材マークの描画*/
    AGR(fp,r1,r2,t1,t2,dr1,dr2);
    fprintf(fp,"#\n");
    /*構造線描画*/
    switch(iii){
        case 0:
            fprintf(fp,"plot %.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin); /*鉄管内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin+t0);   /*鉄管外径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle\n",rout);        /*掘削径*/
            fprintf(fp,"#\n");
            SYMB_CRACK(fp,n,R,L,a);/*ひび割れ描画*/
            break;
        case 1:
            fprintf(fp,"plot %.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin); /*水路内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle\n",rout);        /*水路外径*/
            fprintf(fp,"#\n");
            break;
        case 2:
            fprintf(fp,"plot %.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin); /*水路内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin+t0);   /*内側鉄筋外径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle\n",rout);        /*水路外径*/
            fprintf(fp,"#\n");
            break;
        case 3:
            fprintf(fp,"plot %.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin); /*水路内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin+t0);   /*内側鉄筋外径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle\n",rout);        /*水路外径*/
            fprintf(fp,"#\n");
            SYMB_CRACK(fp,n,R,L,a);/*ひび割れ描画*/
            break;
        case 4:
            fprintf(fp,"plot %.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin); /*水路内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle,\\\n",rin+t0);   /*内側鉄筋外径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle,\\\n",rout-t0);  /*外側鉄筋内径*/
            fprintf(fp,"%.3f with lines linewidth 5 linetype 1 notitle\n",rout);        /*水路外径*/
            fprintf(fp,"#\n");
            SYMB_CRACK(fp,n,R,L,a);/*ひび割れ描画*/
            break;
    }
    fprintf(fp,"#\n");
    SYMB_ROCK(fp,rout);/*岩盤記号描画*/
    fprintf(fp,"#\n");
    fprintf(fp,"unset multiplot\n");
    fprintf(fp,"#\n");
    fclose(fp);

}
/*---------------------------------------------------------------------------*/
void STR_PS(FILE *fp,int iii,double rin,double rout,double t0)
{
    /*埋設鉄管：鉄管，コンクリート，岩盤の文字*/
    double x1,y1,x2,y2,dr,phi,rr=2.0;
    char str[50];

    if(iii==0){
    /*鉄管内径，板厚，コンクリート外径用矢印*/
    fprintf(fp,"set style arrow 1 size graph 0.02,20,50 filled linewidth 1\n");
        phi=20.0;
        x1=0.0;
        y1=0.0;
        x2=rin*cos(phi/180.0*pi);
        y2=rin*sin(phi/180.0*pi);
    fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
        phi=20.0;
        x1=(rin+t0+3.0*t0)*cos(phi/180.0*pi);
        y1=(rin+t0+3.0*t0)*sin(phi/180.0*pi);
        x2=(rin+t0)*cos(phi/180.0*pi);
        y2=(rin+t0)*sin(phi/180.0*pi);
    fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
        phi=-30.0;
        x1=0.0;
        y1=0.0;
        x2=rout*cos(phi/180.0*pi);
        y2=rout*sin(phi/180.0*pi);
    fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
    fprintf(fp,"#\n");
    /*鉄管，コンクリート，岩盤の文字*/
    strcpy(str,"Steel pipe");
        x1=0.0;
        y1=-(rin-1.2);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    strcpy(str,"Filling");
        x1=0.0;
        y1=-(rin+t0)-0.25*(rout-rin-t0);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    strcpy(str,"concrete");
        x1=0.0;
        y1=-(rin+t0)-0.6*(rout-rin-t0);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    strcpy(str,"Bed rock");
        x1=0.0;
        y1=-rout-1.5;
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    fprintf(fp,"#\n");
    /*鉄管内径，板厚，コンクリート外径*/
    dr=0.5;
    strcpy(str,"r{/=12 s}");
        phi=20.0;
        x1=0.5*rin*cos(phi/180.0*pi)-dr*sin(phi/180.0*pi);
        y1=0.5*rin*sin(phi/180.0*pi)+dr*cos(phi/180.0*pi);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center rotate by %.0f\n",str,x1,y1,phi);
    strcpy(str,"t");
        phi=20.0;
        x1=(rin+t0+2.0*t0)*cos(phi/180.0*pi)-dr*sin(phi/180.0*pi);
        y1=(rin+t0+2.0*t0)*sin(phi/180.0*pi)+dr*cos(phi/180.0*pi);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center rotate by %.0f\n",str,x1,y1,phi);
    strcpy(str," r{/=12 c}");
        phi=-30.0;
        x1=rout*cos(phi/180.0*pi);
        y1=rout*sin(phi/180.0*pi);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f left rotate by %.0f\n",str,x1,y1,phi);
    fprintf(fp,"#\n");
    }

}
/*---------------------------------------------------------------------------*/
void STR_RC(FILE *fp,int iii,double rin,double rout,double t0)
{
    /*RC水路：コンクリート，岩盤の文字*/
    double x1,y1,x2,y2,phi=70.0,rr=2.0;
    char str[50];

    strcpy(str,"Concrete");
        x1=0.0;
        y1=-rin-0.4*(rout-rin);
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    strcpy(str,"Bed rock");
        x1=0.0;
        y1=-rout-1.5;
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
    fprintf(fp,"#\n");

    /* 鉄筋表示文字 */
    if(2<=iii&&iii<=3){
        strcpy(str,"Single reinforcement");
            x1=(rout+rr)*cos(phi/180.0*pi);
            y1=(rout+rr)*sin(phi/180.0*pi)+0.6;
            fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
        fprintf(fp,"set style arrow 1 size graph 0.02,20,50 filled linewidth 1\n");
        x1=(rout+rr)*cos(phi/180.0*pi);
        y1=(rout+rr)*sin(phi/180.0*pi);
        x2=(rin+t0)*cos(phi/180.0*pi);
        y2=(rin+t0)*sin(phi/180.0*pi);
        fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
        fprintf(fp,"#\n");
    }
    if(iii==4){
        strcpy(str,"Double reinforcements");
            x1=(rout+rr)*cos(phi/180.0*pi);
            y1=(rout+rr)*sin(phi/180.0*pi)+0.6;
            fprintf(fp,"set label \"%s\" at %.3f,%.3f center\n",str,x1,y1);
        fprintf(fp,"set style arrow 1 size graph 0.02,20,50 filled linewidth 1\n");
        x1=(rout+rr)*cos(phi/180.0*pi);
        y1=(rout+rr)*sin(phi/180.0*pi);
        x2=(rout)*cos(phi/180.0*pi);
        y2=(rout)*sin(phi/180.0*pi);
        fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
        x1=(rout)*cos(phi/180.0*pi);
        y1=(rout)*sin(phi/180.0*pi);
        x2=(rin+t0)*cos(phi/180.0*pi);
        y2=(rin+t0)*sin(phi/180.0*pi);
        fprintf(fp,"set arrow as 1 from %.3f,%.3f to %.3f,%.3f\n",x1,y1,x2,y2);
        fprintf(fp,"#\n");
    }
}
/*---------------------------------------------------------------------------*/
void AGR(FILE *fp,double r1,double r2,double t1,double t2,double dr1,double dr2)
{
    /*骨材マークの描画*/
    int i,MG,MT;
    char str[10];
    double x[Nmax],y[Nmax];

    strcpy(str,"o");
    CAL_XYARRAY(r1,r2,t1,t2,dr1,dr2,x,y,&MG,&MT);
    for(i=0;i<=MG;i++){
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center font \"arial,%d\"\n",str,x[i],y[i],15);
    }
    for(i=MG+1;i<=MT;i++){
        fprintf(fp,"set label \"%s\" at %.3f,%.3f center font \"arial,%d\"\n",str,x[i],y[i],5);
    }
}
/*---------------------------------------------------------------------------*/
void SYMB_CRACK(FILE *fp,int n,double R,double L,double a)
{
    /*ひび割れ描画*/
    int i,j,m;
    double x[22],y[22];
    double dL,co_r,co_t,phi;

    fprintf(fp,"set origin 0,0\n");
    fprintf(fp,"plot \"-\" with lines linewidth 5 linetype 1 notitle,\\\n");   /*ひび割れ：0度*/
    fprintf(fp,"\"-\" with lines linewidth 5 linetype 1 notitle,\\\n");        /*ひび割れ：60度*/
    fprintf(fp,"\"-\" with lines linewidth 5 linetype 1 notitle,\\\n");        /*ひび割れ：120度*/
    fprintf(fp,"\"-\" with lines linewidth 5 linetype 1 notitle,\\\n");        /*ひび割れ：180度*/
    fprintf(fp,"\"-\" with lines linewidth 5 linetype 1 notitle,\\\n");        /*ひび割れ：240度*/
    fprintf(fp,"\"-\" with lines linewidth 5 linetype 1 notitle\n");           /*ひび割れ：280度*/
    fprintf(fp,"#\n");

    dL=L/(double)n;
    m=2*(n+1)-1;
    x[0]=R;
    for(i=1;i<=m-1;i++){
        x[i]=R+0.25*dL*(double)(2*i-1);
    }
    x[m]=R+L;
    y[0]=0.0;
    for(i=1;i<=m-2;i=i+2){
        y[i]=a;
    }
    for(i=2;i<=m-1;i=i+2){
        y[i]=-a;
    }
    y[m]=0.0;

    for(i=0;i<=5;i++){
        phi=60.0*(double)i;
        for(j=0;j<=m;j++){
            CAL_TR(phi,R,x[j],y[j],&co_r,&co_t);
            fprintf(fp,"%.3f %.3f\n",co_t,co_r);
            if(j==m)fprintf(fp,"e\n");
        }
    }
}
/*---------------------------------------------------------------------------*/
void SYMB_ROCK(FILE *fp,double rout)
{
    /*岩盤記号描画*/
    int i;
    double theta;
    double d_theta=15.0;

    for(i=0;i<=3;i++){
        switch(i){
            case 0:theta=  0.0;break;
            case 1:theta= 90.0;break;
            case 2:theta=180.0;break;
            case 3:theta=270.0;break;
        }
        fprintf(fp,"set origin 0,0\n");
        fprintf(fp,"set samples 120\n");
        fprintf(fp,"set trange [%.0f:%.0f]\n",theta-d_theta,theta+d_theta);
        fprintf(fp,"plot %.3f-0.5*sqrt(abs(sin(t/5.0*180.0))) with lines linewidth 5 linetype 1 notitle\n",rout+0.5);
    }

}
/*---------------------------------------------------------------------------*/
void CAL_TR(double phi,double R,double x1,double y1,double *co_r,double *co_t)
{
    /****************************/
    /* ひび割れ点の座標変換     */
    /****************************/
    double x2,y2;

    x2=x1*cos(phi/180.0*pi)-y1*sin(phi/180.0*pi);
    y2=x1*sin(phi/180.0*pi)+y1*cos(phi/180.0*pi);
    *co_r=sqrt(x2*x2+y2*y2);
    if(x2==0.0&&0.0<=y2)*co_t=90.0;
    if(x2==0.0&&y2<=0.0)*co_t=270.0;
    if(x2!=0.0)*co_t=atan(y2/x2)/pi*180.0;;
    if(90.0<phi&&phi<=270.0)*co_t=180.0+*co_t;
}
/*---------------------------------------------------------------------------*/
void CAL_XY(double r1,double r2,double t1,double t2,double *xx,double *yy)
{
    /****************************/
    /* 乱数での座標設定         */
    /****************************/
    double rr,tt;

    rr=r1+(double)rand()/(double)RAND_MAX*(r2-r1);
    tt=t1+(double)rand()/(double)RAND_MAX*(t2-t1);
    *xx=rr*cos(tt/180.0*pi);
    *yy=rr*sin(tt/180.0*pi);
}
/*---------------------------------------------------------------------------*/
void CAL_XYARRAY(double r1,double r2,double t1,double t2,double dr1,double dr2,double x[Nmax],double y[Nmax],int *MG,int *MT)
{
    /****************************/
    /* 骨材を描画する座標の設定 */
    /****************************/
    int i,k,l1,l2,mflag,iflag;
    double xx,yy;

    CAL_XY(r1,r2,t1,t2,&xx,&yy);
    x[0]=xx;
    y[0]=yy;
    k=0;

    l1=0;
    do{
        l2=0;
        do{
            CAL_XY(r1,r2,t1,t2,&xx,&yy);
            mflag=1;
            for(i=0;i<=k;i++){
                if(dr1<sqrt((xx-x[i])*(xx-x[i])+(yy-y[i])*(yy-y[i]))){
                    iflag=1;
                }else{
                    iflag=0;
                }
                mflag=mflag*iflag;
            }
            if(mflag==1)break;
            l2=l2+1;
        }while(l2<(int)(Nmax/2));
        if(mflag==1){
            k=k+1;
            x[k]=xx;
            y[k]=yy;
        }
        l1=l1+1;
    }while(l1<(int)(Nmax/2));
    *MG=k;

    l1=0;
    do{
        l2=0;
        do{
            CAL_XY(r1,r2,t1,t2,&xx,&yy);
            mflag=1;
            for(i=0;i<=k;i++){
                if(dr2<sqrt((xx-x[i])*(xx-x[i])+(yy-y[i])*(yy-y[i]))){
                    iflag=1;
                }else{
                    iflag=0;
                }
                mflag=mflag*iflag;
            }
            if(mflag==1)break;
            l2=l2+1;
        }while(l2<(int)(Nmax/2));
        if(mflag==1){
            k=k+1;
            x[k]=xx;
            y[k]=yy;
        }
        l1=l1+1;
    }while(l1<(int)(Nmax/2));
    *MT=k;

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