/*---------------------------*/
/* プログラム：gccHIST.c */
/*---------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define ND 30000
#define NK 100
/*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,
    int ndata,double dmean,double dsig,double dmax,double dmin);
void QSORT(int n,double data[ND]);
void CALC(int n,double data[ND],double xmin,double xmax,double ds,
    int kdata,double plotx[NK],double ploty[NK]);
void MinMax(int n,double data[ND],
    double *dmean,double *dmin,double *dmax,double *dsig);
/*---------------------------------------------------------------------------*/

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,kdata;
    double pdata[ND],ds,dmean,dmin,dmax,dsig,px,alpha;
    double plotx[NK],ploty[NK];
    struct GPdata gp={
        "","","",
        "",0.0,0.0,0.0,"",
        "",0.0,0.0,0.0,""
    };
    FILE *fin,*fout;

    /*データ入力*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncat(gp.scom,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        strcpy(sv,strtok(dat,","));
        strncat(gp.sxjiku,sv,strlen(sv)-1);
    fgets(dat,sizeof dat,fin);
        kdata=atoi(strtok(dat,","))-1;  /*分割数-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,","));
            strncat(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);
            pdata[i]=atof(strtok(dat,","));
    }
    fclose(fin);

    strcpy(gp.syjiku,"Relative frequency (%)");
    strcpy(gp.slogy,"N");
    gp.ymin=0.0;
    gp.ymax=100.0;
    gp.dy=10.0;
    strcpy(gp.syfmt,"\"%.0f\"");

    /*対数軸変換*/
    if(strncmp(gp.slogx,"L",1)==0){
        for(i=0;i<=ndata;i++){
            if(pdata[i]<gp.xmin)pdata[i]=gp.xmin;
            pdata[i]=log10(pdata[i]);
        }
        gp.xmin=log10(gp.xmin);
        gp.xmax=log10(gp.xmax);
    }

    ds=(gp.xmax-gp.xmin)/(double)(kdata+1);
    QSORT(ndata,pdata);
    CALC(ndata,pdata,gp.xmin,gp.xmax,ds,kdata,plotx,ploty);
    MinMax(ndata,pdata,&dmean,&dmin,&dmax,&dsig);

    /*対数軸の普通軸変換*/
    if(strncmp(gp.slogx,"L",1)==0){
        for(i=0;i<=ndata;i++){
            pdata[i]=pow(10.0,pdata[i]);
        }
        gp.xmin=pow(10.0,gp.xmin);
        gp.xmax=pow(10.0,gp.xmax);
        for(i=0;i<=kdata;i++){
            plotx[i]=pow(10.0,plotx[i]);
        }
        dmean=pow(10.0,dmean);
        dmin=pow(10.0,dmin);
        dmax=pow(10.0,dmax);
        dsig=pow(10.0,dsig);
    }

    /*ヒストグラム一時データファイル書き出し*/
    fout=fopen("temp0.prn","w");
    for(i=0;i<=kdata;i++){
        fprintf(fout,"%g %g\n",plotx[i],ploty[i]/(double)(ndata+1)*100.0);
    }
    fclose(fout);

    /*累積確率一時データファイル書き出し*/
    alpha=0.5; /*Hazen plot*/
    fout=fopen("temp1.prn","w");
    for(i=0;i<=ndata;i++){
        j=i+1;
        px=(j-alpha)/((double)(ndata+1)+1.0-2.0*alpha)*100.0;
        fprintf(fout,"%g %g\n",pdata[i],px);
    }
    fclose(fout);

    /*gnuplotによる作図*/
    GPLOT(fnameW,gp,ndata,dmean,dsig,dmax,dmin);

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

void GPLOT(char fnameW[50],struct GPdata gp,
    int ndata,double dmean,double dsig,double dmax,double dmin)
{
    FILE *fp;
    char sv[20]="",dat[50]="";

    /****************************/
    /*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);
        /*ｘ軸定義*/
        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);
        }
        /*ｙ軸定義*/
        if(strncmp(gp.slogy,"L",1)==0)fprintf(fp,"set logscale y\n");
        if(strncmp(gp.slogy,"N",1)==0){
            fprintf(fp,"set ytics %f\n",gp.dy);
            fprintf(fp,"set format y %s\n",gp.syfmt);
        }
        /*軸範囲定義*/
        fprintf(fp,"set xrange [%f:%f]\n",gp.xmin,gp.xmax);
        fprintf(fp,"set yrange [%f:%f]\n",gp.ymin,gp.ymax);
        /*ラベル定義*/
        sprintf(sv,"%d",ndata+1);strcpy(dat,"n=");strcat(dat,sv);
        fprintf(fp,"set label \"%s\" at graph 0.02,graph 0.95 left\n",dat);
        sprintf(sv,"%g",dmean);strcpy(dat,"mean=");strcat(dat,sv);
        fprintf(fp,"set label \"%s\" at graph 0.02,graph 0.90 left\n",dat);
        sprintf(sv,"%g",dsig);strcpy(dat,"sigma=");strcat(dat,sv);
        fprintf(fp,"set label \"%s\" at graph 0.02,graph 0.85 left\n",dat);
        sprintf(sv,"%g",dmax);strcpy(dat,"max=");strcat(dat,sv);
        fprintf(fp,"set label \"%s\" at graph 0.02,graph 0.80 left\n",dat);
        sprintf(sv,"%g",dmin);strcpy(dat,"min=");strcat(dat,sv);
        fprintf(fp,"set label \"%s\" at graph 0.02,graph 0.75 left\n",dat);
        /*軸グリッド描画設定*/
        fprintf(fp,"set grid xtics ytics mxtics mytics\n");
        /*塗りつぶし設定*/
        fprintf(fp,"set style fill solid 0.2 border -1\n");
        /*プロット実行（プロットは各種定義が終わって最後に実行）*/
        fprintf(fp,"plot \\\n");
        fprintf(fp,"\"temp0.prn\" with boxes linewidth 1 linetype 5 notitle, \\\n");
        fprintf(fp,"\"temp1.prn\" with lines linewidth 1 linetype 1 notitle\n");
    fclose(fp);

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

void QSORT(int n,double data[ND])
{
    /* n：データ個数-1 */
    int nn,k,sw,i,j;
    double di,dj;
    nn=n+1;
    k=nn;
CONT1:
    k=(int)(k/2);
CONT2:
    sw=0;
    for(i=0;i<=nn-k-1;i++){
        j=i+k;
        di=data[i];dj=data[j];
        if(data[i]>data[j]){
            data[i]=dj;data[j]=di;sw=1;
        }
    }
    if(sw==1)goto CONT2;
    if(k>1)goto CONT1;

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

void CALC(int n,double data[ND],double xmin,double xmax,double ds,int kdata,double plotx[NK],double ploty[NK])
{
int k,num,i;
double xx1,xx2;

    /* n：データ個数-1 */
    xx1=xmin;
    for(k=0;k<=kdata;k++){
        xx2=xx1+ds;
        num=0;
        plotx[k]=xx1+0.5*ds;
        for(i=0;i<=n-1;i++){
            if(xx1<=data[i]&&data[i]<xx2)num=num+1;
        }
        ploty[k]=(double)num;
        xx1=xx2;
    }
}
/*---------------------------------------------------------------------------*/

void MinMax(int n,double data[ND],double *dmean,double *dmin,double *dmax,double *dsig)
{
    int i;
    /* n：データ個数-1 */
    *dmean=0.0;
    for(i=0;i<=n;i++){
        *dmean=*dmean+data[i];
    }
    *dmean=*dmean/(double)(n+1);
    *dsig=0;
    for(i=0;i<=n;i++){
        *dsig=*dsig+(data[i]-*dmean)*(data[i]-*dmean);
    }
    *dsig=sqrt(*dsig/(double)n);
    *dmin=data[0];
    *dmax=data[n];
}
/*---------------------------------------------------------------------------*/