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

my <- mean(vec_lq)
sy <- sd(vec_lq)
ry <- nn*sum((vec_lq-my)^3)/(nn-1)/(nn-2)/sy^3

bb<-4/ry/ry
aa<-sy/sqrt(bb)
cc<-my-aa*bb
vec_x <- exp(cc+aa*qgamma(vec_p,shape=bb,rate=1.0))
rr <- cor(vec_x,vec_q)
#
ss1 <- (log(vec_q)-cc)/aa
ss2 <- qgamma(vec_p,shape=bb,rate=1.0)
ww  <- sum((ss1-ss2)^2)
s99 <- qgamma(0.99,shape=bb,rate=1.0)
s01 <- qgamma(0.01,shape=bb,rate=1.0)
SLSC<- sqrt(ww/nn)/abs(s99-s01)
#
postscript(file="fig_R_HF_qq_LP3_mom.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","LP3 (Moment)")
s1 <- sprintf(" n=%g\n",nn)
s2 <- sprintf(" r=%.3f\n",rr)
s3 <- sprintf(" a=%.3f\n",aa)
s4 <- sprintf(" b=%.3f\n",bb)
s5 <- sprintf(" c=%.3f\n",cc)
s6 <- sprintf(" SLSC=%.3f\n",SLSC)
mtext(text=s0,side=3,line=-3,adj=0)
mtext(text=s1,side=3,line=-4.5,adj=0)
mtext(text=s2,side=3,line=-6,adj=0)
mtext(text=s3,side=3,line=-7.5,adj=0)
mtext(text=s4,side=3,line=-9,adj=0)
mtext(text=s5,side=3,line=-10.5,adj=0)
mtext(text=s6,side=3,line=-12.0,adj=0)
dev.off()