提交: 66d8f97d0eb7dffd173df61cf2ea11e93828d0e1 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: 66d8f97d0eb7dffd173df61cf2ea11e93828d0e1 (tree)
• 时间: 2008-04-28 19:36:04
• 作者: iselllo
• Commiter: iselllo

Log Message

I added the code inset.R which shows how to create and inset plot using
R.

更改概述

• add: R-codes/inset.R

差异

diff -r 9933ec258969 -r 66d8f97d0eb7 R-codes/inset.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/R-codes/inset.R	Mon Apr 28 10:36:04 2008 +0000
@@ -0,0 +1,187 @@
+rm(list=ls())
+
+
+#x<-0:10;
+time <- read.table("time_red.dat", header=FALSE)
+time <- as.matrix(time)
+#y<-x^4;
+delta_r<- read.table("delta_r_sq_time.dat", header=FALSE)
+delta_r <- as.matrix(delta_r)
+
+
+#now a linear fit
+
+my_fit <- lm(delta_r[4:length(delta_r)]~ time[4:length(time)])
+
+b <- as.real(summary(my_fit)$coefficients[ ,1][2])
+
+print ("the slope is, ")
+
+print (b)
+
+a <- as.real(summary(my_fit)$coefficients[ ,1][1])
+
+print ("the intercept is, ")
+print (a)
+
+fit1 <- a+b*time[4:length(time)]
+
+#print ("my_fit is, ")
+#summary(my_fit)
+
+time_short<- read.table("time_red_short.dat", header=FALSE)
+time_short<- as.matrix(time_short)
+#y<-x^4;
+delta_r_short<- read.table("delta_r_sq_time_short.dat", header=FALSE)
+delta_r_short<- as.matrix(delta_r_short)
+
+
+
+time_inset <- c(time_short[2:length(time_short)],time[5:15])
+
+delta_r_inset <- c(delta_r_short[2:length(delta_r_short)], delta_r[5:15])
+
+
+
+
+my_fit_short<- lm(log10(delta_r_short[2:10])~  log10(time_short[2:10]))
+
+
+b_short<- as.real(summary(my_fit_short)$coefficients[ ,1][2])
+
+print ("the slope is, ")
+
+print (b_short)
+
+a_short<- as.real(summary(my_fit_short)$coefficients[ ,1][1])
+
+print ("the intercept is, ")
+print (a_short)
+
+
+#fit_short <- a_short+b_short*log10(time_short[2:length(time_short)])
+
+
+fit_short <- a_short+b_short*log10(time_inset)
+
+
+
+
+my_fit_log<- lm(log10(delta_r[4:length(delta_r)])~ log10(time[4:length(time)]))
+
+b <- as.real(summary(my_fit_log)$coefficients[ ,1][2])
+
+print ("the slope is, ")
+
+print (b)
+
+a <- as.real(summary(my_fit_log)$coefficients[ ,1][1])
+
+print ("the intercept is, ")
+print (a)
+
+#fit_log<- a+b*log10(time_short[2:length(time_short)])
+
+fit_log<- a+b*log10(time_inset)
+
+
+
+library(gridBase)
+#X11(width=8,height=8)
+
+# produce outer ("main") plot
+pdf("inset_plot.pdf")
+   par( mar = c(4.5,5, 2, 1) + 0.1)
+
+   plot(time,delta_r,col="blue",lwd=1.5,lty=1,pch=1,xlim=range(c(0,400)),
+        ylim=range(c(0,25)),
+        xlab=expression("Time"),ylab=expression("<"* delta*""*r ^2 *">"),cex.axis=1.4,cex.lab=1.6)
+   lines(time[4:length(time)],fit1,col="black",lwd=2.)
+
+legend(350,10,cex=1.2, c(expression("Simulation"),
+
+                   expression("Linear fit"),
+               expression("Power-law fit")), 
+lwd=c(2,2,2), lty=c(-1,1,2),pch = c(1,-1,-1),col=c("blue","black","red"),box.lwd=0,box.lty=0,
+,xjust = 1, yjust = 1)
+
+
+
+vp <- baseViewports()
+   pushViewport(vp$inner,vp$figure,vp$plot)
+
+# push viewport that will contain the inset
+
+#   pushViewport(viewport(x=0.1,y=0.9,width=.5,height=.5,just=c("left","top")))
+
+   pushViewport(viewport(x=-0.11,y=1.06,width=.68,height=.65,just=c(0,1)))
+
+
+# now either define viewport to contain the whole inset figure
+
+#   par(fig=gridFIG(),new=T) # or gridPLT() # ...or just the plotting are (coordinate system)
+
+   par(fig=gridPLT(),new=T)
+
+#   par(plt=gridPLT(),new=T)
+
+
+# draw frame around selected area (for illustration only)
+
+#   grid.rect(gp=gpar(lwd=3,col="red"))
+
+# plot inset figure
+
+##    plot(time_short,delta_r_short,"b",col="blue",xlab="",ylab="",ylim=range(c(0,max(delta_r_short+0.02)))
+##         ,yaxt="n",cex.axis=1.4,cex.lab=1.6)
+##    axis(side=2, at=c( 0, 0.075, 0.15),
+##    labels=expression(0, 0.075, 0.15),cex.axis=1.4,cex.lab=1.6)
+
+##    plot(time_short[2:length(time_short)],delta_r_short[2:length(delta_r_short)],
+##         pch=1,col="blue",xlab="",ylab=""
+##         ,cex.axis=1.4,cex.lab=1.6,log="xy")
+
+   plot(time_inset,delta_r_inset,
+        pch=1,col="blue",xlab="",ylab=""
+        ,cex.axis=1.4,cex.lab=1.6,log="xy",xaxt="n",yaxt="n")
+    axis(side=1, at=c( 0.1,1, 5,15),
+    labels=expression(0.1,1, 5,15),cex.axis=1.4,cex.lab=1.6)
+
+    axis(side=2, at=c( 0.00002,0.01,0.5),
+    labels=expression(2%*%10^{-5},10^{-1},0.5),cex.axis=1.4,cex.lab=1.6)
+
+
+
+##    axis(side=2, at=c( log10(0.01), log10(0.075), log10(0.15)),
+##    labels=expression(log10(0.01), log10(0.075), log10(0.15)),cex.axis=1.4,cex.lab=1.6)
+
+
+##    lines(time_short[2:length(time_short)],10**(fit_short), col="black", lwd=2.)
+
+##    lines(time_short[2:length(time_short)],10**(fit_log), col="red", lwd=2.)
+
+   lines(time_inset,10**(fit_short), col="red", lwd=2.,lty=2)
+
+   lines(time_inset,10**(fit_log), col="black", lwd=2.)
+
+
+
+   #axis(side=2, at=c( 0, 0.075, 0.15),
+   #labels=expression(0, 0.075, 0.15),cex.axis=1.4,cex.lab=1.6)
+
+
+
+
+
+# pop all viewports from stack
+popViewport(4)
+
+dev.off()
+
+
+pdf("log-log_plot")
+   plot(time_short[2:length(time_short)],delta_r_short[2:length(delta_r_short)],"b",col="blue"
+        ,cex.axis=1.4,cex.lab=1.6,log="xy")
+dev.off()
+
+print ("So far so good")