/*----------------------------*/
/* プログラム：gccHAZEN-png.c */
/*----------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define ND 1000
/*gnuplot用構造体宣言*/
/*scom,sxjiku,syjiku :タイトル(凡例),x軸ラベル,ｙ軸ラベル*/
/*slogx,xmin,xmax,dx,sxfmt :x軸種,x軸最小値,x軸最大値,x軸増分,x軸書式*/
/*slogy,ymin,ymax,dy,syfmt :y軸種,y軸最小値,y軸最大値,y軸増分,y軸書式*/
struct GPdata{
    char scom[100];char sxjiku[100];char syjiku[100];
    char slogx[3];double xmin;double xmax;double dx;char sxfmt[20];
    char slogy[3];double ymin;double ymax;double dy;char syfmt[20];
};
/*関数プロトタイプ*/
void main_part(char fnameR[50],char fnameW[50]);
void GPLOT(char fnameW[50],struct GPdata gp,double aa,double bb,double rr);
double TODA(double rx);
void KAIKI(int ndata, double datax[ND], double datay[ND],
    double *aa, double *bb, double *rr);
/*---------------------------------------------------------------------------*/

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[50],char fnameW[50])
{
    char dat[256]="",sv[100]="";
    int i,j,ndata;
    double datax[ND],datay[ND],work;
    double aa,bb,rr;
    struct GPdata gp={
        "","","",
        "",0.0,0.0,0.0,"",
        "",0.0,0.0,0.0,""
    };
    FILE *fin,*fout;
    double alpha = 0.5; /*Hazen公式*/

    /*データ入力*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncpy(gp.scom,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncpy(gp.sxjiku,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        strcpy(gp.slogx,strtok(dat,","));
        gp.xmin=atof(strtok(NULL,","));
        gp.xmax=atof(strtok(NULL,","));
        if(strncmp(gp.slogx,"N",1)==0){
            gp.dx=atof(strtok(NULL,","));
            strcpy(sv,strtok(NULL,","));
            strncpy(gp.sxfmt,sv,strlen(sv)-1);
        }
    fgets(dat,sizeof dat,fin);
        ndata=atoi(strtok(dat,","))-1;
    for(i=0;i<=ndata;i++){
        fgets(dat,sizeof dat,fin);
            datax[i]=atof(strtok(dat,","));
    }
    fclose(fin);

    /*ｘ軸値ソート(大きい順)*/
    for(i=0;i<=ndata;i++){
        for(j=i+1;j<=ndata;j++){
            if(datax[i]<datax[j]){
                work=datax[i];
                datax[i]=datax[j];
                datax[j]=work;
            }
        }
    }

    /*ｘ軸対数の場合の常用対数変換*/
    if(strncmp(gp.slogx,"L",1)==0){
        for(i=0;i<=ndata;i++){
            if(datax[i]<gp.xmin)datax[i]=gp.xmin;
            datax[i]=log10(datax[i]);
        }
    }

    /*超過確率計算とｙ軸データ処理*/
    for(i=0;i<=ndata;i++){
        datay[i]=((double)(i+1)-alpha)/((double)(ndata+1)+1.0-2.0*alpha);
        work=datay[i];
        datay[i]=TODA(work);
    }

    /*回帰係数算定（y=aa*x+bb [相関係数：rr]）*/
    KAIKI(ndata,datax,datay,&aa,&bb,&rr);

    /*ｘ軸対数値を普通軸数値に再変換*/
    if(strncmp(gp.slogx,"L",1)==0){
        for(i=0;i<=ndata;i++){
            datax[i]=pow(10.0,datax[i]);
        }
    }

    /*一時データファイル書き出し*/
    fout=fopen("temp.prn","w");
    for(i=0;i<=ndata;i++){
        fprintf(fout,"%f %f\n",datax[i],datay[i]);
    }
    fclose(fout);

    /*gnuplotによる作図*/
    GPLOT(fnameW,gp,aa,bb,rr);

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

void GPLOT(char fnameW[50],struct GPdata gp,double aa,double bb,double rr)
{
    FILE *fp;
    int i;
    char str[20]="";
    char seq[100]="";
    char dat[100]="";
    double yy,work;

    strcpy(gp.syjiku,"Exceedance probability (%)");
    /****************************/
    /*gnuplot制御用ファイル作成*/
    /****************************/
    fp=fopen("gpl_temp.txt","w");
        /*出力ファイル定義*/
        fprintf(fp,"set terminal png\n");
        fprintf(fp,"set output \"%s\"\n",fnameW);
        /*グラフタイトル定義*/
        fprintf(fp,"set title \"%s\"\n",gp.scom);
        /*軸ラベル定義*/
        fprintf(fp,"set xlabel \"%s\"\n",gp.sxjiku);
        fprintf(fp,"set ylabel \"%s\"\n",gp.syjiku);
        /*y軸（確率軸）定義*/
        fprintf(fp,"set ytics ( \\\n");
        strcpy(str,"99.9");work=atof(str)/100.0;gp.ymax=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,gp.ymax);
        strcpy(str,"99")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"95")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"90")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"80")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"50")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"20")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"10")  ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"5")   ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"1")   ;work=atof(str)/100.0;yy=TODA(work);fprintf(fp,"\"%s\" %f,\\\n",str,yy);
        strcpy(str,"0.1") ;work=atof(str)/100.0;gp.ymin=TODA(work);fprintf(fp,"\"%s\" %f)\n",str,gp.ymin);
        /*ｘ軸定義*/
        if(strncmp(gp.slogx,"L",1)==0){
            fprintf(fp,"set logscale x\n");
            fprintf(fp,"set mxtics 10\n");
        }
        if(strncmp(gp.slogx,"N",1)==0){
            fprintf(fp,"set xtics %f\n",gp.dx);
            fprintf(fp,"set format x %s\n",gp.sxfmt);
        }
        /*軸範囲定義*/
        fprintf(fp,"set xrange [%f:%f]\n",gp.xmin,gp.xmax);
        fprintf(fp,"set yrange [%f:%f]\n",gp.ymin,gp.ymax);
        /*回帰直線設定*/
        fprintf(fp,"a=%g\n",aa);
        fprintf(fp,"b=%g\n",bb);
        if(strncmp(gp.slogx,"L",1)==0)strcpy(seq,"f(x)=a*log10(x)+b");
        if(strncmp(gp.slogx,"N",1)==0)strcpy(seq,"f(x)=a*x+b");
        fprintf(fp,"%s\n",seq);
        /*回帰係数ラベル設定*/
        sprintf(str,"a=%g",aa);strcpy(dat,str);strcat(dat,", ");
        sprintf(str,"b=%g",bb);strcat(dat,str);strcat(dat,", ");
        sprintf(str,"r=%.3f",rr);strcat(dat,str);
        fprintf(fp,"set label \"%s\" at graph %g,graph %g\n",dat,0.0,1.03);
        /*軸グリッド描画設定*/
        fprintf(fp,"set grid xtics ytics mxtics mytics\n");
        /*プロット実行（プロットは各種定義が終わって最後に実行）*/
        fprintf(fp,"plot f(x) with lines title \"%s\", \\\n",seq);
        fprintf(fp,"\"temp.prn\" with points pointtype 4 linetype 4 notitle\n");
    fclose(fp);

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

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

double TODA(double rx)
{
    /* 非超過確率の％点 */
    double ry,ay,uay;
    double b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10;

    if(rx==0.5){
        uay=0.0;
    }else{
        if(rx<0.5){
            ry=rx;
        }else{
            ry=1-rx;
        }
        ay=-log(4.0*ry*(1.0-ry));
        b0 = 1.570796288;
        b1 = 0.03706987906;
        b2 =-0.0008364353589;
        b3 =-0.0002250947176;
        b4 = 0.000006841218299;
        b5 = 0.000005824238515;
        b6 =-0.00000104527497;
        b7 = 0.00000008360937017;
        b8 =-0.000000003231081277;
        b9 = 0.00000000003657763036;
        b10= 0.0000000000006936233982;
        uay=sqrt(ay*(b0+b1*ay+b2*pow(ay,2.0)+b3*pow(ay,3.0)+b4*pow(ay,4.0)+b5*pow(ay,5.0)+b6*pow(ay,6.0)+b7*pow(ay,7.0)+b8*pow(ay,8.0)+b9*pow(ay,9.0)+b10*pow(ay,10.0)));
        if(rx<0.5)uay=-uay;
    }
    return uay;
}
/*---------------------------------------------------------------------------*/

void KAIKI(int ndata, double datax[ND], double datay[ND],
           double *aa, double *bb, double *rr)
{
    int i;
    double x1,y1,x2,xy,xm,ym,c1,c2,c3;
            
    /* 回帰式：y=aa*x+bb */
    x1 = 0.0 ; y1 = 0.0 ; x2 = 0.0 ; xy = 0.0;
    for(i=0;i<=ndata;i++){
        x1 = x1 + datax[i];
        y1 = y1 + datay[i];
        x2 = x2 + datax[i] * datax[i];
        xy = xy + datax[i] * datay[i];
    }
    xm = x1 / (double)(ndata + 1) ; ym = y1 / (double)(ndata + 1);
    *aa = ((double)(ndata + 1) * xy - x1 * y1) / ((double)(ndata + 1) * x2 - x1 * x1);
    *bb = (x2 * y1 - x1 * xy) / ((double)(ndata + 1) * x2 - x1 * x1);
    c1 = 0.0 ; c2 = 0.0 ; c3 = 0.0;
    for(i=0;i<=ndata;i++){
        c1 = c1 + (datax[i] - xm) * (datay[i] - ym);
        c2 = c2 + (datax[i] - xm) * (datax[i] - xm);
        c3 = c3 + (datay[i] - ym) * (datay[i] - ym);
    }
    *rr = c1 / sqrt(c2 * c3);
}
/*---------------------------------------------------------------------------*/