/*----------------------------------------------------------------------------*/
/* gccSCVFsurge.c */
/*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define g 9.8    /* 重力加速度 */
#define pi M_PI  /* 円周率 */
#define N 10     /* 水槽断面検討ケース数 */

char title[100];        /* タイトル            （入力値）*/
double Hg;              /* 総落差              （入力値）*/
double Q0;              /* 最大使用水量        （入力値）*/
double L;               /* 圧力水路延長        （入力値）*/
double c;               /* 圧力水路損失水頭係数（入力値）*/
double Cd;              /* 制水口流量係数      （入力値）*/
double zm;              /* 最大上昇水位 (規制値：入力値) */
double d0;              /* 圧力水路内径        （入力値）*/
double xc;              /* 選定したポート径 */
double yc;              /* 選定した立坑径 */
double zc;              /* 選定した立坑径の最高水位 */
double R_DFC2;          /* 動的安定性より定まる最小水槽径(記憶用) */
double DF[N];           /* 水槽内径（計算条件格納） */
double zmax[N];         /* 最大上昇水位（計算結果格納） */

double draw_xmin,draw_xmax,draw_dx;
double draw_ymin,draw_ymax,draw_dy;

void INP_DAT(char fnameR[50]);
void CAL_SC();
void CAL_VF();
double CAL_ZM(double d1,double k0,double h0,double md,int *iflag);
void PLOT_SC(char fnameW_SC[50]);
void PLOT_VF(char fnameW_VF[50]);
/*----------------------------------------------------------------------------*/

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

    strcpy(fnameR,argv[1]);
    strcpy(fnameW_SC,argv[2]);
    strcpy(fnameW_VF,argv[3]);

    INP_DAT(fnameR);
    CAL_SC();
    PLOT_SC(fnameW_SC);
    CAL_VF();
    PLOT_VF(fnameW_VF);

    printf("End of calculation.");

    return 0;
}
/*----------------------------------------------------------------------------*/
void INP_DAT(char fnameR[])
{
    FILE *fp;
    char dat[100];

    fp=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fp);
        strncpy(title,dat,strlen(dat)-1); /*タイトル*/
    fgets(dat,sizeof dat,fp);
    fgets(dat,sizeof dat,fp);
        Hg =atof(strtok(dat,","));        /* 総落差 */
        Q0 =atof(strtok(NULL,","));       /* 最大使用水量 */
        L  =atof(strtok(NULL,","));       /* 圧力水路延長 */
        d0 =atof(strtok(NULL,","));       /* 圧力水路内径 */
        c  =atof(strtok(NULL,","));       /* 圧力水路損失水頭係数 */
        Cd =atof(strtok(NULL,","));       /* 制水口流量係数 */
        zm =atof(strtok(NULL,","));       /* 最大上昇水位 */
        xc =atof(strtok(NULL,","));       /* 選定ポート径 */
        yc =atof(strtok(NULL,","));       /* 選定立坑径 */
    fclose(fp);

}
/*----------------------------------------------------------------------------*/
void CAL_SC()
{
    FILE *fp;
    int iflag;
    double fwa,v0;
    double h0;                 /* 遮断前圧力水路損失水頭 h0=c*v0*v0 */
    double k0;                 /* 遮断時制水口抵抗 */
    double DFC1,DFC2,DFC3;     /* 水槽内径(限界条件) */
    double Fsa;                /* 水槽断面積 */
    double Fpa;                /* 制水口断面積 */
    double d2,d2min,d2max,d2d; /* 制水口内径計算範囲指定 */

    /* 安定条件 */
    fwa=0.25*pi*d0*d0;
    v0=Q0/fwa;
    h0=c*v0*v0;

    Fsa=L*fwa/c/2.0/g/(Hg-zm);
    DFC2=sqrt(4.0*Fsa/pi);
    R_DFC2=DFC2;

    d2min=0.1;
    d2max=(double)((int)(DFC2+1.0));
    d2d=0.1;
 
    fp=fopen("temp_SC.prn","w");
    printf("d2  k0-h0  DFC1  DFC2  DFC3\n");
    fprintf(fp,"#d2  k0-h0  DFC1  DFC2  DFC3\n");
    d2=d2min-d2d;
    iflag=-1;
    do{
        d2=d2+d2d;
        Fpa=0.25*pi*d2*d2;
        k0=0.5/g*Q0*Q0/Cd/Cd/Fpa/Fpa;
        if(k0-h0<0)break;
 
        Fsa=h0*L*fwa/c/(h0+k0)/g/Hg;
        DFC1=sqrt(4.0*Fsa/pi);
 
        Fsa=h0*L*fwa*fwa*fwa/2.0/g/c/c/k0/Q0/Q0;
        DFC3=sqrt(4.0*Fsa/pi);
 
        printf("%6.3f  %10.3e  %6.3f  %6.3f  %6.3f\n",d2,k0-h0,DFC1,DFC2,DFC3);
        fprintf(fp,"%6.3f  %10.3e  %6.3f  %6.3f  %6.3f\n",d2,k0-h0,DFC1,DFC2,DFC3);
    }while(d2<d2max);
    fclose(fp);

    draw_xmin=0.0;
    draw_xmax=(double)((int)(d2-d2d));
    draw_dx=0.5;
    draw_ymin=0.0;
    draw_ymax=(double)((int)((DFC2+1.0)/2.0))*6.0;
    draw_dy=2.0;
}
/*----------------------------------------------------------------------------*/
void CAL_VF()
{
    FILE *fp;
    int i,iflag;
    char dat[256];
    char sv[20];
    double fwa,v0,md;
    double Fsa;                /* 水槽断面積 */
    double Fpa;                /* 制水口内径,断面積 */
    double d1,d1min,d1d;       /* 水槽内径計算範囲指定 */
    double d2,d2min,d2max,d2d; /* 制水口内径計算範囲指定 */
    double h0;                 /* 遮断前圧力水路損失水頭 h0=c*v0*v0 */
    double k0;                 /* 遮断時制水口抵抗 */

    /* Vogt-Forchheimer の式による最高上昇水位 */
    fwa=0.25*pi*d0*d0;
    v0=Q0/fwa;
    h0=c*v0*v0;

    d1min=(double)((int)(R_DFC2+1.0));
    d1d=1.0;
    for(i=0;i<=N-1;i++){
        DF[i]=d1min+d1d*(double)i;
    }
    fp=fopen("temp_VF.prn","w");
    printf("===========================================================================\n");
    strcpy(dat,"#d2 k0-h0");
    for(i=0;i<=N-1;i++){
        strcat(dat,"  DF=");sprintf(sv,"%.1f",DF[i]);strcat(dat,sv);
    }
    printf("%s\n",dat);
    fprintf(fp,"%s\n",dat);
    d2min=0.1;
    d2d=0.1;
    d2=d2min-d2d;
    do{
        d2=d2+d2d;
        Fpa=0.25*pi*d2*d2;
        k0=0.5/g*Q0*Q0/Cd/Cd/Fpa/Fpa;
        if(k0-h0<0)break;
        iflag=1; /* log関数の引数を監視するフラグ */
        for(i=0;i<=N-1;i++){
            d1=DF[i];
            Fsa=0.25*pi*d1*d1;
            md=2.0*g*Fsa*(h0+k0)/L/fwa/v0/v0;
            zmax[i]=CAL_ZM(d1,k0,h0,md,&iflag);
        }
        if(iflag==1){
            sprintf(sv,"%8.3f",d2);strcpy(dat,sv);
            sprintf(sv,"%8.3f",k0-h0);strcat(dat," ");strcat(dat,sv);
            for(i=0;i<=N-1;i++){
                sprintf(sv,"%8.3f",zmax[i]);strcat(dat," ");strcat(dat,sv);
            }
            printf("%s\n",dat);
            fprintf(fp,"%s\n",dat);
        }
    }while(d2<draw_xmax);
    fclose(fp);

    iflag=1;
    d1=yc;  /*水槽内径*/
    d2=xc;  /*制水口内径*/
    Fsa=0.25*pi*d1*d1;
    Fpa=0.25*pi*d2*d2;
    k0=0.5/g*Q0*Q0/Cd/Cd/Fpa/Fpa;
    md=2.0*g*Fsa*(h0+k0)/L/fwa/v0/v0;
    zc=CAL_ZM(d1,k0,h0,md,&iflag); /*選択諸元での最高上昇水位*/

    fp=fopen("temp_VF.prn","r");
    fgets(dat,sizeof dat,fp);
    while(fgets(dat,sizeof dat,fp)!=NULL){
        strtok(dat," ");
        strtok(NULL," ");
        for(i=0;i<=N-1;i++){
            zmax[i]=atof(strtok(NULL," "));
        }
    }
    fclose(fp);
    draw_ymin=0.0;
    draw_ymax=(double)((int)(0.1*zmax[0])+1)*10.0;
    draw_dy=10.0;

    draw_ymax=100.0;

}
/*----------------------------------------------------------------------------*/
double CAL_ZM(double d1,double k0,double h0,double md,int *iflag)
{
    double ff,df,z1,z2,z;
    double eps=1e-10;       /* 収束計算0判定 */

    if(md*k0<1.0){
        z1=-1.0/md+0.0001;
        do{
            z=z1;
            if(1.0+md*z<0){
                *iflag=0;break;
            }
            ff=((1.0+md*z)-log(1.0+md*z))-((1.0+md*h0)-log(1.0-md*k0));
            df=md*(1.0-1.0/(1.0+md*z));
            z2=z1-ff/df;
            z1=z2;
        }while(eps<fabs(ff));
    }
    if(1.0<md*k0){
        z1=1.0/md+0.0001;
        do{
            z=z1;
            ff=((md*z-1.0)+log(md*z-1.0))-(log(md*k0-1.0)-(md*h0+1.0));
            if(md*z-1.0<0){
                *iflag=0;break;
            }
            df=md*(1.0+1.0/(md*z-1.0));
            z2=z1-ff/df;
            z1=z2;
        }while(eps<fabs(ff));
    }

    return fabs(z1);
}
/*----------------------------------------------------------------------------*/
void PLOT_SC(char fnameW_SC[50])
{
    FILE *fp;
    int i;
    double xx,yy;

    /****************************/
    /*gnuplot制御用ファイル作成*/
    /****************************/
    fp=fopen("gpl_temp_SC.txt","w");
    fprintf(fp,"set terminal postscript eps enhanced font \"Helvetica\" 16\n");
    fprintf(fp,"set output \"%s\"\n",fnameW_SC);
    fprintf(fp,"set size ratio 1.0\n");
    if(0<strstr(title,"Load"))fprintf(fp,"set key left top\n");
    if(0<strstr(title,"Input"))fprintf(fp,"set key right top\n");
    fprintf(fp,"set key width -2\n");
    fprintf(fp,"set key spacing 1.2\n");
    fprintf(fp,"set key box\n");
    fprintf(fp,"set xlabel \"Diameter of the port (m)\"\n");
    fprintf(fp,"set ylabel \"Diameter of the surge tank (m)\"\n");
    fprintf(fp,"set xrange [%.3f:%.3f]\n",draw_xmin,draw_xmax);
    fprintf(fp,"set xtics %.3f\n",draw_dx);
    fprintf(fp,"set format x \"%s\"\n","%.1f");
    fprintf(fp,"set yrange [%.3f:%.3f]\n",draw_ymin,draw_ymax);
    fprintf(fp,"set ytics %.3f\n",draw_dy);
    fprintf(fp,"set format y \"%s\"\n","%.0f");
    fprintf(fp,"set tmargin 2\n");
    fprintf(fp,"set label \"  %s\" at graph 0.0,graph 1.03 left\n",title);
    fprintf(fp,"set rmargin 8\n");
    /*2本の直線(nohead矢印)による選択点の表示*/
    fprintf(fp,"set style arrow 2 size graph 0.05,20 filled nohead linewidth 1 linetype 1\n");
    fprintf(fp,"set arrow as 2 from %.3f,%.3f to %.3f,%.3f\n",draw_xmin,yc,xc,yc);
    fprintf(fp,"set arrow as 2 from %.3f,%.3f to %.3f,%.3f\n",xc,draw_ymin,xc,yc);
    xx=(xc-draw_xmin)/(draw_xmax-draw_xmin)+0.02;
    yy=0.0;
    fprintf(fp,"set label \"  %.1fm\" at graph %.3f,graph %.3f left rotate by 90\n",xc,xx,yy);
    xx=0.0;
    yy=(yc-draw_ymin)/(draw_ymax-draw_ymin)+0.02;
    fprintf(fp,"set label \"  %.1fm\" at graph %.3f,graph %.3f left\n",yc,xx,yy);
    /*動的安定性を満足する水槽径最小値の表示*/
    xx=0.0;
    yy=(R_DFC2-draw_ymin)/(draw_ymax-draw_ymin)+0.02;
    fprintf(fp,"set label \"   %.3fm (Target Z{/=10 m}=%.1fm)\" at graph %.3f,graph %.3f left\n",R_DFC2,zm,xx,yy);
    /*軸作成*/
    fprintf(fp,"set grid xtics ytics mxtics mytics\n");
    fprintf(fp,"plot \"temp_SC.prn\" using 1:3 with lines linewidth 3 linetype 4 title \"D{/=10 req} > D{/=10 1}\", \\\n");
    fprintf(fp,"     \"temp_SC.prn\" using 1:4 with lines linewidth 3 linetype 1 title \"D{/=10 req} > D{/=10 2}\", \\\n");
    fprintf(fp,"     \"temp_SC.prn\" using 1:5 with lines linewidth 3 linetype 2 title \"D{/=10 req} < D{/=10 3}\", \\\n");
    fprintf(fp,"     \"-\" with points pointsize 2 pointtype 7 title \"Selected point\"\n");
    fprintf(fp,"     %.3f %.3f\n",xc,yc);
    fprintf(fp,"     e\n");
    fclose(fp);

    /****************************/
    /*gnuplot制御*/
    /****************************/
    fp=popen("gnuplot.exe","w");
        fprintf(fp,"load \"gpl_temp_SC.txt\"\n");
        fflush(fp);
    pclose(fp);

}
/*----------------------------------------------------------------------------*/
void PLOT_VF(char fnameW_VF[50])
{
    FILE *fp;
    int i;
    double xx,yy;

    /****************************/
    /*gnuplot制御用ファイル作成*/
    /****************************/
    fp=fopen("gpl_temp_VF.txt","w");
    fprintf(fp,"set terminal postscript eps enhanced font \"Helvetica\" 16\n");
    fprintf(fp,"set output \"%s\"\n",fnameW_VF);
    fprintf(fp,"set size ratio 1.0\n");
    if(0<strstr(title,"Load"))fprintf(fp,"set key left top\n");
    if(0<strstr(title,"Input"))fprintf(fp,"set key right top\n");
    fprintf(fp,"set key width -4\n");
    fprintf(fp,"set key spacing 1.2\n");
    fprintf(fp,"set key box\n");
    fprintf(fp,"set xlabel \"Diameter of the port (m)\"\n");
    fprintf(fp,"set ylabel \"Maximum WL  Z{/=10 max} (m)  or  Difference of pressure  |Z'{/=10 rm}| (m)\"\n");
    fprintf(fp,"set xrange [%.3f:%.3f]\n",draw_xmin,draw_xmax);
    fprintf(fp,"set xtics %.3f\n",draw_dx);
    fprintf(fp,"set format x \"%s\"\n","%.1f");
    fprintf(fp,"set yrange [%.3f:%.3f]\n",draw_ymin,draw_ymax);
    fprintf(fp,"set ytics %.3f\n",draw_dy);
    fprintf(fp,"set format y \"%s\"\n","%.0f");
    fprintf(fp,"set tmargin 2\n");
    fprintf(fp,"set label \"  %s\" at graph 0.0,graph 1.03 left\n",title);
    fprintf(fp,"set rmargin 8\n");
    /*直線(nohead矢印)による水位上昇限界値の表示*/
    fprintf(fp,"set style arrow 2 size graph 0.05,20 filled nohead linewidth 1 linetype 1\n"); 
    fprintf(fp,"set arrow as 2 from %.3f,%.3f to %.3f,%.3f\n",draw_xmin,zm,draw_xmax,zm);
    xx=0.0;
    yy=(zm-draw_ymin)/(draw_ymax-draw_ymin)+0.02;
    fprintf(fp,"set label \"  Target Z{/=10 m}=%.1fm\" at graph %.3f,graph %.3f left\n",zm,xx,yy);
    /*最高上昇水位内径凡例*/
    for(i=0;i<=N-1;i++){
        fprintf(fp,"set label \" D=%.0fm\"  at first %.3f,first %.3f left\n",DF[i],draw_xmax,zmax[i]);
    }
    /*軸作成*/
    fprintf(fp,"set grid xtics ytics mxtics mytics\n");
    fprintf(fp,"plot \\\n");
    for(i=0;i<=N-1;i++){
        fprintf(fp,"\"temp_VF.prn\" using 1:%d with line linewidth 3 linetype 1 notitle, \\\n",3+i);
    }
    fprintf(fp,"\"temp_VF.prn\" using 1:2 with linespoints linewidth 3 pointtype 6 title \"|Z'{/=10 rm}| = k{/=10 0} - h{/=10 0}\", \\\n");
    fprintf(fp,"\"-\" with points pointsize 2 pointtype 7 title \"Selected point\"\n");
    fprintf(fp,"%.3f %.3f\n",xc,zc);
    fprintf(fp,"e\n");
    fclose(fp);
    /****************************/
    /*gnuplot制御*/
    /****************************/
    fp=popen("gnuplot.exe","w");
        fprintf(fp,"load \"gpl_temp_VF.txt\"\n");
        fflush(fp);
    pclose(fp);

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