/*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define NDMAX 1000

void main_part(char fnameR[],char fnameW[]);
void TeX_DRN_XJIKU(FILE *fout,double xmin,double xmax,double dx,
    double kxi,double kxf,double kyi,double kyf);
void TeX_DRL_XJIKU(FILE *fout,double xmin,double xmax,
    double kxi,double kxf,double kyi,double kyf);
void TeX_DRH_YJIKU(FILE *fout,double ymin,double ymax,
    double kxi,double kxf,double kyi,double kyf);
double TODA(double rx);
void KAIKI(int nd,double datax[NDMAX],double datay[NDMAX],
    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[],char fnameW[])
{
    /*作図用変数*/
    double HSIZE;/*画像横寸法（mm）*/
    double DHL;/*画像左余白寸法：縦軸描画用（mm）*/
    double DHR;/*画像右余白寸法（mm）*/
    double VSIZE;/*画像縦寸法（mm）*/
    double DVL;/*画像下余白寸法：横軸描画用（mm）*/
    double DVU;/*画像上余白寸法（mm）*/
    char sxjiku[50]="";/*ｘ軸名（単位含む）*/
    double xmin;/*ｘ軸最小値*/
    double xmax;/*ｘ軸最大値*/
    double dx;/*ｘ軸増分*/
    char slogx[2]="";/*ｘ軸の種類（L：対数，N：普通）*/
    char syjiku[50]="";/*ｙ軸名*/
    double ymin;/*ｙ軸最小値*/
    double ymax;/*ｙ軸最大値*/
    double xx,yy,kxx,kyy,kxx1,kyy1,kxx2,kyy2,dd;
    double kxi,kxf,kyi,kyf;
    int nd;/*データ数−１*/
    double alpha=0.5;/*Hazen公式*/
    double aa,bb,rr;/*回帰係数および相関係数*/
    /*プロットデータ格納配列宣言*/
    double datax[NDMAX];
    double datay[NDMAX];
    int i,j,k;
    double work;
    char strs[50];
    /*入出力用変数宣言*/
    FILE *fin,*fout;
    char dat[256];
    char sv[10];

    /*データ読み込み*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
        strcpy(sxjiku,strtok(dat,","));
    fgets(dat,sizeof dat,fin);
        HSIZE=atof(strtok(dat,","));
        DHL  =atof(strtok(NULL,","));
        DHR  =atof(strtok(NULL,","));
    fgets(dat,sizeof dat,fin);
        VSIZE=atof(strtok(dat,","));
        DVL  =atof(strtok(NULL,","));
        DVU  =atof(strtok(NULL,","));
    fgets(dat,sizeof dat,fin);
        xmin=atof(strtok(dat,","));
        xmax=atof(strtok(NULL,","));
        dx  =atof(strtok(NULL,","));
        strcpy(slogx,strtok(NULL,","));
    fgets(dat,sizeof dat,fin);
        nd=atof(strtok(dat,","))-1;
    for(i=0;i<=nd;i++){
        fgets(dat,sizeof dat,fin);
            datax[i]=atof(strtok(dat,","));
    }
    fclose(fin);

    /*プロット領域定義*/
    kxi=DHL;
    kxf=HSIZE-DHR;
    kyi=DVL;
    kyf=VSIZE-DVU;

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

    /*ｘ軸対数値の場合0データ処理*/
    if(strncmp(slogx,"L",1)==0){
        for(i=0;i<=nd;i++){
            if(datax[i]<pow(10.0,xmin)){datax[i]=pow(10.0,xmin);}
            datax[i]=log10(datax[i]);
        }
    }

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

    /*回帰式(y=aa*x+bb,rr)*/
    KAIKI(nd,datax,datay,&aa,&bb,&rr);

    /*TeXファイル書き込み*/
    fout=fopen(fnameW,"w");
    fprintf(fout,"\\documentclass[a4paper,10pt]{jsarticle}\n");
    fprintf(fout,"\\input{jartfmt}\n");
    fprintf(fout,"\\usepackage[dvipdfm]{graphics,graphicx}\n");
    fprintf(fout,"\\usepackage{colortbl}\n");
    fprintf(fout,"\\begin{document}\n");
    fprintf(fout,"\\begin{center}\n");
    fprintf(fout,"\\unitlength=1mm\n");
    fprintf(fout,"\\begin{picture}(%.3f,%.3f)\n",HSIZE,VSIZE);
    /*ｘ軸描画*/
    if(strncmp(slogx,"N",1)==0){TeX_DRN_XJIKU(fout,xmin,xmax,dx,kxi,kxf,kyi,kyf);}
    if(strncmp(slogx,"L",1)==0){TeX_DRL_XJIKU(fout,xmin,xmax,kxi,kxf,kyi,kyf);}
    /*ｙ軸描画*/
    ymax=TODA(0.999);
    ymin=TODA(0.001);
    TeX_DRH_YJIKU(fout,ymin,ymax,kxi,kxf,kyi,kyf);
    /*枠線描画*/
    fprintf(fout,"\\path(%.3f,%.3f)(%.3f,%.3f)\n",kxi,kyi,kxf,kyi);
    fprintf(fout,"\\path(%.3f,%.3f)(%.3f,%.3f)\n",kxf,kyi,kxf,kyf);
    fprintf(fout,"\\path(%.3f,%.3f)(%.3f,%.3f)\n",kxf,kyf,kxi,kyf);
    fprintf(fout,"\\path(%.3f,%.3f)(%.3f,%.3f)\n",kxi,kyf,kxi,kyi);
    /*ｘ軸名描画*/
    fprintf(fout,"\\put(%.3f,%.3f){\\makebox(0,0){%s}}\n",0.5*(kxi+kxf),kyi-8.0,sxjiku);
    /*ｙ軸名描画*/
    strcpy(syjiku,"超過確率~(\\%)");
    fprintf(fout,"\\put(%.3f,%.3f){\\makebox(0,0){\\rotatebox{90}{%s}}}\n",kxi-15.0,0.5*(kyi+kyf),syjiku);
    /*丸印描画*/
    for(i=0;i<=nd;i++){
        xx=datax[i];
        yy=datay[i];
        if(yy<ymin)yy=ymin;
        kxx=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
        kyy=kyi+(yy-ymin)*(kyf-kyi)/(ymax-ymin);
        strcpy(strs,"$\\bullet$");
        fprintf(fout,"\\put(%.3f,%.3f){\\makebox(0,0){%s}}\n",kxx,kyy,strs);
    }
    /*回帰曲線描画*/
    xx=xmin;yy=aa*xx+bb;
    if(yy>=ymax){
    yy=ymax;xx=(yy-bb)/aa;
    }
    kxx1=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
    kyy1=kyi+(yy-ymin)*(kyf-kyi)/(ymax-ymin);
    xx=xmax;yy=aa*xx+bb;
    if(yy<ymin){
    yy=ymin;xx=(yy-bb)/aa;
    }
    kxx2=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
    kyy2=kyi+(yy-ymin)*(kyf-kyi)/(ymax-ymin);
    fprintf(fout,"\\path(%.3f,%.3f)(%.3f,%.3f)\n",kxx1,kyy1,kxx2,kyy2);
    /*回帰係数・相関係数表示*/
    dd=2.0;
    fprintf(fout,"\\put(%.3f,%.3f){\\begin{minipage}{2.5cm}\n",kxf-25.0,kyf-dd*5.0);
    fprintf(fout,"\\footnotesize\n");
    fprintf(fout,"\\begin{tabular}{|l|}\\hline\n");
    fprintf(fout,"\\rowcolor{white} $y=a\\cdot x+b$\\\\ \n");
    fprintf(fout,"\\rowcolor{white} $a=%.5f$\\\\ \n",aa);
    fprintf(fout,"\\rowcolor{white} $b=%.5f$\\\\ \n",bb);
    fprintf(fout,"\\rowcolor{white} $r=%.5f$\\\\ \\hline \n",rr);
    fprintf(fout,"\\end{tabular}\n");
    fprintf(fout,"\\end{minipage}}\n");

    fprintf(fout,"\\end{picture}\n");
    fprintf(fout,"\\end{center}\n");
    fprintf(fout,"\\end{document}\n");

    fclose(fout);

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

void TeX_DRN_XJIKU(FILE *fout,double xmin,double xmax,double dx,
    double kxi,double kxf,double kyi,double kyf)
{
    int i,ix;
    double xx,kxx,wv;
    char sv[10];
    /*普通ｘ軸描画*/
    ix=(int)((xmax-xmin)/dx);
    for(i=0;i<=ix;i++){
        xx=xmin+(double)i*dx;
        kxx=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
        wv=xmin+(double)i*dx;
        sprintf(sv,"%.0f",wv);
        if((int)(dx*1000.0)%10!=0){sprintf(sv,"%.3f",wv);}
        if((int)(dx*1000.0)%10==0&&(int)(dx*100.0)%10!=0){sprintf(sv,"%.2f",wv);}
        if((int)(dx*1000.0)%10==0&&(int)(dx*100.0)%10==0&&(int)(dx*10.0)%10!=0){sprintf(sv,"%.1f",wv);}
        fprintf(fout,"\\dottedline{1}(%.3f,%.3f)(%.3f,%.3f)\n",kxx,kyi,kxx,kyf);
        fprintf(fout,"\\put(%.3f,%.3f){\\makebox(0,0){%s}}\n",kxx,kyi-3.0,sv);
    }
}
/*---------------------------------------------------------------------------------*/

void TeX_DRL_XJIKU(FILE *fout,double xmin,double xmax,
    double kxi,double kxf,double kyi,double kyf)
{
    int i,j;
    double xx,kxx;
    /*対数ｘ軸描画*/
    for(i=(int)xmin;i<=(int)xmax-1;i++){
        for(j=1;j<=9;j++){
            xx=log10((double)j*pow(10.0,(double)i));
            kxx=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
            fprintf(fout,"\\dottedline{1}(%.3f,%.3f)(%.3f,%.3f)\n",kxx,kyi,kxx,kyf);
        }
    }
    for(i=(int)xmin;i<=(int)xmax;i++){
        xx=(double)i;
        kxx=kxi+(xx-xmin)*(kxf-kxi)/(xmax-xmin);
        fprintf(fout,"\\put(%.3f,%.3f){\\makebox(0,0){10$^{%d}$}}\n",kxx,kyi-3.0,i);
    }
}
/*---------------------------------------------------------------------------------*/

void TeX_DRH_YJIKU(FILE *fout,double ymin,double ymax,
    double kxi,double kxf,double kyi,double kyf)
{
    int i;
    double yy,kyy,work;
    char strs[10];

    for(i=1;i<=11;i++){
        if(i== 1)strcpy(strs,"99.9");
        if(i== 2)strcpy(strs,"99");
        if(i== 3)strcpy(strs,"95");
        if(i== 4)strcpy(strs,"90");
        if(i== 5)strcpy(strs,"80");
        if(i== 6)strcpy(strs,"50");
        if(i== 7)strcpy(strs,"20");
        if(i== 8)strcpy(strs,"10");
        if(i== 9)strcpy(strs,"5");
        if(i==10)strcpy(strs,"1");
        if(i==11)strcpy(strs,"0.1");
        work=(double)(atof(strs)/100.0);
        yy=TODA(work);
        kyy=kyi+(yy-ymin)*(kyf-kyi)/(ymax-ymin);
        fprintf(fout,"\\dottedline{0.5}(%.3f,%.3f)(%.3f,%.3f)\\put(%.3f,%.3f){\\makebox(0,0){%s}}\n",kxi,kyy,kxf,kyy,kxi-5.0,kyy,strs);
    }
}
/*---------------------------------------------------------------------------------*/

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 nd,double datax[NDMAX],double datay[NDMAX],
    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<=nd;i++){
        x1=x1+datax[i];
        y1=y1+datay[i];
        x2=x2+datax[i]*datax[i];
        xy=xy+datax[i]*datay[i];
    }
    xm=x1/(double)(nd+1);ym=y1/(double)(nd+1);
    *aa=((double)(nd+1)*xy-x1*y1)/((double)(nd+1)*x2-x1*x1);
    *bb=(x2*y1-x1*xy)/((double)(nd+1)*x2-x1*x1);
    c1=0.0;c2=0.0;c3=0.0;
    for(i=0;i<=nd;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);
}
/*---------------------------------------------------------------------------------*/