d <- read.table('inp_RF_M.txt',skip=1) 
vec_q <- sort(d[,1])
nn <- length(vec_q)
alp <- 0.4 #Cunnane formula
vec_p <- ppoints(nn,alp)
#plot(vec_q,vec_p)

aa <- (vec_q[1]*vec_q[nn]-median(vec_q)^2)/(vec_q[1]+vec_q[nn]-2*median(vec_q))
vec_y <- log(vec_q-aa)
my <- mean(vec_y)
sy <- sd(vec_y)*sqrt((nn-1)/nn)
vec_x <- aa+exp(my+sy*qnorm(vec_p))
rr <- cor(vec_x,vec_q)
#
ss1 <- (log(vec_q-aa)-my)/sy
ss2 <- qnorm(vec_p)
ww  <- sum((ss1-ss2)^2)
s99 <- qnorm(0.99)
s01 <- qnorm(0.01)
SLSC<- sqrt(ww/nn)/abs(s99-s01)
#
postscript(file="fig_R_HF_qq_LN3_IWAI.eps", width=10,family="Helvetica")
par(ps=18)
par(pty="s")
plot(vec_q,vec_x,log="xy",xlim=c(20,500),ylim=c(20,500),xlab="Observed Q",ylab="Estimated Q",pch=1,cex=1.2)
abline(0,1)
s0 <- sprintf(" %s\n","LN3 (IWAI)")
s1 <- sprintf(" n=%g\n",nn)
s2 <- sprintf(" r=%.3f\n",rr)
s3 <- sprintf(" a=%.3f\n",aa)
s4 <- sprintf(" my=%.3f\n",my)
s5 <- sprintf(" sy=%.3f\n",sy)
s6 <- sprintf(" SLSC=%.3f\n",SLSC)
mtext(text=s0,side=3,line=-3.0,adj=0)
mtext(text=s1,side=3,line=-4.5,adj=0)
mtext(text=s2,side=3,line=-6.0,adj=0)
mtext(text=s3,side=3,line=-7.5,adj=0)
mtext(text=s4,side=3,line=-9.0,adj=0)
mtext(text=s5,side=3,line=-10.5,adj=0)
mtext(text=s6,side=3,line=-12.0,adj=0)
dev.off()