#**********************************
# Logarithmic normal distribution
#**********************************
{
    if(2<=NR)datax[NR-1]=$1
}
END{
    nd=NR-1
    #--------------------------
    # Sort (small->large)
    #--------------------------
    for(i=1;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;
            }
        }
    }
    #--------------------------
    # Non-exceeding probability
    # Plotting and positioned formula
    #    Weibull : alpha=0
    #    Cunnane : alpha=0.4
    #    Hazen   : alpha=0.5
    #--------------------------
    alpha=0.5
    for(i=1;i<=nd;i++){
        p=(i-alpha)/(nd+1-2*alpha)
        datay[i]=TODA(p,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10)
    }
    #--------------------------
    # Regression (y=aa*x+bb)
    # (logarithmic normal distribution)
    #--------------------------
    for(i=1;i<=nd;i++){
        datax[i]=log(datax[i])/log(10)
    }
    x1=0.0;y1=0.0;x2=0.0;xy=0.0;
    for(i=1;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/nd;ym=y1/nd;
    aa=(nd*xy-x1*y1)/(nd*x2-x1*x1);
    bb=(x2*y1-x1*xy)/(nd*x2-x1*x1);
    c1=0.0;c2=0.0;c3=0.0;
    for(i=1;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);
    #--------------------------
    # Print out
    #--------------------------
    print  "#y=a*x+b"
    printf "#a=%4.3e\n",aa
    printf "#b=%4.3e\n",bb
    printf "#r=%4.3f\n",rr
    printf "#nd=%d\n",nd
    print  "#"
    print  "#"
    print  "#*** Regression line ***"
    print  "#",10^datax[1],aa*datax[1]+bb
    print  "#",10^datax[nd],aa*datax[nd]+bb
    print  "#*** y-axis data (text) ***"
    nep=0.999;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.995;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.99 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.95 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.90 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.50 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.10 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.05 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.01 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.005;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    nep=0.001;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#",-0.02,yy,14,0,0,"MR",nep
    print "#*** y-axis data (grid) ***"
    nep=0.999;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.995;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.99 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.95 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.90 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.50 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.10 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.05 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.01 ;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.005;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    nep=0.001;yy=TODA(nep,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10);print "#>";printf "#0 %g\n",yy;printf "#1 %g\n",yy
    print "#*** Measured data ***"
    for(i=1;i<=nd;i++)printf "%g %g\n",10^datax[i],datay[i]
}
function TODA(rx,     ry,ay,uay,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10)
{
    # --------------------
    # Toda's approxmation
    # --------------------
    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*ay^2+b3*ay^3+b4*ay^4+b5*ay^5+b6*ay^6+b7*ay^7+b8*ay^8+b9*ay^9+b10*ay^10));
        if(rx<0.5)uay=-uay;
    }
    return uay;
}