#**********************************
# Gumbel 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]=-log(-log(p))
    }
    #--------------------------
    # Regression (y=aa*x+bb)
    #--------------------------
    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);
    theta=1.0/aa
    mu=-bb/aa;
    #--------------------------
    # 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
    printf "#theta=%4.3e\n",theta
    printf "#mu=%4.3e\n",mu
    print  "#*** Regression line ***"
    print  "#",datax[1],aa*datax[1]+bb
    print  "#",datax[nd],aa*datax[nd]+bb
    print  "#*** y-axis data (text) ***"
    p=0.999;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.995;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.99 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.95 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.90 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.80 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.70 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.50 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.10 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.01 ;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    p=0.001;print "#",-0.02,-log(-log(p)),14,0,0,"MR",p
    print "#*** y-axis data (grid) ***"
    p=0.999;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.995;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.99 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.95 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.90 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.80 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.70 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.50 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.10 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.01 ;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    p=0.001;print "#>";printf "#0 %g\n",-log(-log(p));printf "#1 %g\n",-log(-log(p))
    print "#*** Measured data ***"
    for(i=1;i<=nd;i++)printf "%g %g\n",datax[i],datay[i]
}