Commit: 04dd78bb33cee72e8b260c65262ca2d42e46c73c - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 04dd78bb33cee72e8b260c65262ca2d42e46c73c (tree)
• Zeit: 2006-12-11 00:25:45
• Autor: iselllo
• Commiter: iselllo

Log Message

(empty log message)

Ändern Zusammenfassung

• add: codes-for-fitting7/integration-and-plots-corrected-final.R

Diff

diff -r 326dca2a8478 -r 04dd78bb33ce codes-for-fitting7/integration-and-plots-corrected-final.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/codes-for-fitting7/integration-and-plots-corrected-final.R	Sun Dec 10 15:25:45 2006 +0000
@@ -0,0 +1,548 @@
+rm(list=ls())
+
+
+
+
+  #    integral = 0.0D0
+
+
+#      do I=1,NN-4,4
+#!C      S1=h/3.0D0*(f(I)+4.0D0*f(I+2)+f(I+4))
+#!C      S2=h/6.0D0*(f(I)+4.0D0*f(I+1)+2.0D0*f(I+2)+4.0D0*f(I+3)+f(I+4))
+
+
+#!C      S1=h/3.0D0*(f(I)+4.0D0*f(I+1)+2.0D0*f(I+2)+4.0D0*f(I+3)+f(I+4))
+ #                 S2=h/45.0D0*(14.0D0*f(I)+64.0D0*f(I+1)  &
+ #                    &   +24.0D0*f(I+2)+64.0D0*f(I+3)+14.0D0*f(I+4))
+
+
+ #                 integral = integral + S2
+ 
+
+#!C      eps = eps + dabs(S1-S2)/15.0D0
+
+
+#!C      eps = eps + dabs(S1-S2)
+
+
+ #    end do
+
+
+#!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+ #     end subroutine integration
+
+skip<-1
+
+def_integral<-function(integrand,h)
+{
+n<-length(integrand)
+myseq<-seq(1,n-4,4)
+integral<-0
+for (i in myseq)
+{
+integral<- integral + h/45.0*(14.0*integrand[i]+64.0*integrand[i+1] +24.0*integrand[i+2]+64.0*integrand[i+3]+14.0*integrand[i+4])
+
+}
+return(integral)
+}
+
+#Now a test of the integration routine
+
+test_seq<-seq(0,2*pi,len=200)
+integ<-sin(test_seq)
+h2<-test_seq[25]-test_seq[24]
+sin_int<-def_integral(integ,h2)
+
+
+
+lendat<-23
+for(i in 1:lendat)
+{
+ fn <- paste("r",i,".txt",sep="")
+ dat <- read.table(fn,header=FALSE)
+ # ... do your stuff on "dat" here ...
+ dat_temp<-as.matrix(dat)
+{
+if (i ==1) {mymat<-dat_temp}
+
+else
+mymat<-cbind(mymat,dat_temp)
+}
+}
+
+print("OK reading the first files")
+
+mymat_save<-mymat # It is handy to save these data into a separate object and work with mymat only in the following
+
+
+h<-abs(mymat[200,1]-mymat[201,1])
+#vecdim<-dim(mymat)
+#n<-vecdim[2]
+#myseq<-seq(1,n-4,4)
+
+#define below each contribution to the integral and then work out the total.
+
+
+
+
+
+vecdim<-dim(mymat)
+n_int<-vecdim[2]/2
+int_value<-seq(1:n_int)
+int_value2<-seq(1:n_int)
+int_value3<-seq(1:n_int)
+int_value4<-seq(1:n_int)
+int_value5<-seq(1:n_int)
+int_value_save<-seq(1:n_int)
+
+
+
+
+x_seq<-seq(1:n_int)
+
+for (i in 1:n_int) # I define the grid position I will be using
+{
+{ 
+if (i<=11)
+{
+x_seq[i]<-0.01*(i-1)
+}
+
+if (i > 11 & i<=20 )
+{
+x_seq[i]<-0.1*(i-10)
+}
+if (i>20 & i<= lendat)
+{
+x_seq[i]<-i-19
+}
+}
+}
+
+
+for(i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]
+int_value[i]<-def_integral(myseq2,h)
+}
+
+
+
+
+
+print("the density at wall is")
+print(int_value[1])
+
+dens_v<-int_value[ ]/(0.193242)
+
+pdf("density-absorbing-1.pdf")
+plot(x_seq,int_value/int_value[1],col="blue",xlab=expression("x"),ylab=expression(" average v "))
+lines(x_seq,int_value/int_value[1],lwd=2,col="red")
+dev.off()
+
+
+
+save_data<-matrix(ncol=2,nrow=n_int)
+save_data[ ,2]<-int_value[ ]
+save_data[ ,1]<-x_seq[ ]
+
+
+
+
+
+
+write.table(save_data,"density-normalized.txt",quote=FALSE,row.names=FALSE)
+
+#Now I can plot and work out other quantities.
+#mymat[ ,2]<-mymat[ ,1]*mymat[ ,2]
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*mymat[ ,2*(i-1)+1]
+int_value2[i]<-def_integral(myseq2,h)
+}
+
+
+
+pdf("flux-wrong-absorbing-1.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression("J= average (rho*v)"))
+lines(x_seq,int_value2,lwd=2,col="red")
+dev.off()
+
+flux_v<-int_value2
+
+
+
+int_value2[ ]<-int_value2[ ]/int_value[]
+
+pdf("velocity-absorbing-1.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression(" average v"))
+lines(x_seq,int_value2,lwd=2,col="red")
+dev.off()
+
+
+vel_v<-int_value2
+
+flux_v_correct<-vel_v*dens_v
+
+
+
+
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*(mymat[ ,2*(i-1)+1])*(mymat[ ,2*(i-1)+1])
+int_value3[i]<-def_integral(myseq2,h)
+}
+
+
+pdf("rho-velocity-squared-absorbing-1.pdf")
+plot(x_seq,int_value3,col="blue",xlab=expression("x"),ylab=expression(" average [rho*(v squared)]"))
+lines(x_seq,int_value3,lwd=2,col="red")
+dev.off()
+
+
+int_value3[ ]<-int_value3[ ]/int_value[ ] # now int_value3 contains <v>
+
+
+
+pdf("velocity-squared-absorbing-1.pdf")
+plot(x_seq,int_value3,col="blue",xlab=expression("x"),ylab=expression(" average v squared"))
+lines(x_seq,int_value3,lwd=2,col="red")
+dev.off()
+
+velsq_v<-int_value3
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*(mymat[ ,2*(i-1)+1])*(mymat[ ,2*(i-1)+1])
+int_value4[i]<-def_integral(myseq2,h)
+}
+
+pdf("rho-velocity-fluctuations-squared-absorbing-1.pdf")
+plot(x_seq,int_value4,col="blue",xlab=expression("x"),ylab=expression(" average (rho*sig_vv squared)"))
+lines(x_seq,int_value4,lwd=2,col="red")
+dev.off()
+
+
+int_value4[ ]<-int_value4[ ]/int_value[ ]-int_value2[ ]*int_value2[ ] 
+
+
+pdf("velocity-fluctuations-squared-absorbing-1.pdf")
+plot(x_seq,int_value4,col="blue",xlab=expression("x"),ylab=expression(" average (sig_vv squared)"))
+lines(x_seq,int_value4,lwd=2,col="red")
+dev.off()
+
+sig_v<-int_value4
+
+# Now I compare these results to the ones predicted by Alichenkov
+vel_ali1<-function(chi,sig)
+{
+(1-chi)/(1+chi)*sqrt(2*sig/pi)
+}
+
+vel_ali2<-function(chi,vsq)
+{
+(1-chi)/(1+chi)*sqrt(2*vsq/pi)
+}
+
+chi<-0 # perfect absorption
+
+vel_model1<-function()
+{
+vel_ali1(chi=chi,sig=int_value4)
+} 
+
+vel_model2<-function()
+{
+vel_ali2(chi=chi,vsq=int_value3)
+} 
+
+vel1<-vel_model1()
+vel2<-vel_model2()
+
+#Now I add Nagvi's results for the density
+nagvi_den<-function(b,D,x_ini,x)
+{
+b/((D+b*x_ini)*x_ini*(1-0.5*b*x_ini/(D+b*x_ini)))*(x-x_ini)+1/(x_ini*(1-0.5*b*x_ini/(D+b*x_ini)))
+}
+
+nagvi<-function()
+{
+nagvi_den(b,D,x_ini,x)
+}
+
+b<-1
+D<-1/100
+x_ini<-8
+x<-x_seq
+
+ana_dens<-nagvi()
+
+
+
+# Now I'll create a plot with the 3 velocities together
+
+pdf("velocity-absorbing-comparison-Alichenko.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression(" average v"),type="b",
+ylim=range(c(min(int_value2,vel1,vel2),max(int_value2,vel1,vel2))))
+#lines(x_seq,int_value2,lwd=2,col="blue")
+#lines(x_seq,vel1,lwd=2,col="red")
+lines(x_seq,vel1,lwd=2,col="red",type="b")
+#lines(x_seq,vel2,lwd=2,col="black")
+lines(x_seq,vel2,lwd=2,col="black",type="b")
+legend(0.4, 1, c(expression("Numerics"), expression("Alipchenkov eq 6"), expression("Alipchenkov eq 15")), 
+lwd= c(2, 2,2),col=c("blue","red","black" ))
+
+dev.off()
+
+
+###################################################
+##################################################
+
+
+
+rm(mymat)
+rm(dat_temp)
+
+
+for(i in 1:lendat)
+{
+ fn <- paste("../../../no-convection/v0-0/final/rnc",i,".txt",sep="")
+ dat <- read.table(fn,header=FALSE)
+ # ... do your stuff on "dat" here ...
+ dat_temp<-as.matrix(dat)
+{
+if (i ==1) {mymat<-dat_temp}
+
+else
+mymat<-cbind(mymat,dat_temp)
+}
+}
+
+for(i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]
+int_value[i]<-def_integral(myseq2,h)
+}
+
+
+
+
+
+print("the density at wall is")
+print(int_value[1])
+
+dens<-int_value/(0.706045)
+
+pdf("2density-absorbing-1.pdf")
+plot(x_seq,int_value/int_value[1],col="blue",xlab=expression("x"),ylab=expression(" average v "))
+lines(x_seq,int_value/int_value[1],lwd=2,col="red")
+dev.off()
+
+
+
+save_data<-matrix(ncol=2,nrow=n_int)
+save_data[ ,2]<-int_value[ ]
+save_data[ ,1]<-x_seq[ ]
+
+
+
+
+write.table(save_data,"density-normalized.txt",quote=FALSE,row.names=FALSE)
+
+#Now I can plot and work out other quantities.
+#mymat[ ,2]<-mymat[ ,1]*mymat[ ,2]
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*mymat[ ,2*(i-1)+1]
+int_value2[i]<-def_integral(myseq2,h)
+}
+
+
+
+pdf("2flux-absorbing-1.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression("J= average (rho*v)"))
+lines(x_seq,int_value2,lwd=2,col="red")
+dev.off()
+
+flux<-int_value2
+
+int_value2[ ]<-int_value2[ ]/int_value[]
+
+pdf("2velocity-absorbing-1.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression(" average v"))
+lines(x_seq,int_value2,lwd=2,col="red")
+dev.off()
+
+
+vel<-int_value2
+flux_correct<-vel*dens
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*(mymat[ ,2*(i-1)+1])*(mymat[ ,2*(i-1)+1])
+int_value3[i]<-def_integral(myseq2,h)
+}
+
+
+pdf("2rho-velocity-squared-absorbing-1.pdf")
+plot(x_seq,int_value3,col="blue",xlab=expression("x"),ylab=expression(" average [rho*(v squared)]"))
+lines(x_seq,int_value3,lwd=2,col="red")
+dev.off()
+
+
+int_value3[ ]<-int_value3[ ]/int_value[ ] # now int_value3 contains <v>
+
+
+
+pdf("2velocity-squared-absorbing-1.pdf")
+plot(x_seq,int_value3,col="blue",xlab=expression("x"),ylab=expression(" average v squared"))
+lines(x_seq,int_value3,lwd=2,col="red")
+dev.off()
+
+velsq<-int_value3
+
+for (i in 1:n_int)
+{
+myseq2<-mymat[ ,2*i]*(mymat[ ,2*(i-1)+1])*(mymat[ ,2*(i-1)+1])
+int_value4[i]<-def_integral(myseq2,h)
+}
+
+pdf("2rho-velocity-fluctuations-squared-absorbing-1.pdf")
+plot(x_seq,int_value4,col="blue",xlab=expression("x"),ylab=expression(" average (rho*sig_vv squared)"))
+lines(x_seq,int_value4,lwd=2,col="red")
+dev.off()
+
+
+int_value4[ ]<-int_value4[ ]/int_value[ ]-int_value2[ ]*int_value2[ ]
+
+
+pdf("2velocity-fluctuations-squared-absorbing-1.pdf")
+plot(x_seq,int_value4,col="blue",xlab=expression("x"),ylab=expression(" average (sig_vv squared)"))
+lines(x_seq,int_value4,lwd=2,col="red")
+dev.off()
+
+sig<-int_value4
+
+# Now I compare these results to the ones predicted by Alichenkov
+#vel_ali1<-function(chi,sig)
+#{
+#(1-chi)/(1+chi)*sqrt(2*sig/pi)
+#}
+
+#vel_ali2<-function(chi,vsq)
+#{
+#(1-chi)/(1+chi)*sqrt(2*vsq/pi)
+#}
+
+#chi<-0 # perfect absorption
+
+#vel_model1<-function()
+#{
+#vel_ali1(chi=chi,sig=int_value4)
+#} 
+
+#vel_model2<-function()
+#{
+#vel_ali2(chi=chi,vsq=int_value3)
+#} 
+
+vel1bis<-vel_model1()
+vel2bis<-vel_model2()
+
+# Now I'll create a plot with the 3 velocities together
+
+pdf("2velocity-absorbing-comparison-Alichenko.pdf")
+plot(x_seq,int_value2,col="blue",xlab=expression("x"),ylab=expression(" average v"),type="b",
+ylim=range(c(min(int_value2,vel1bis,vel2bis),max(int_value2,vel1bis,vel2bis))))
+#lines(x_seq,int_value2,lwd=2,col="blue")
+#lines(x_seq,vel1,lwd=2,col="red")
+lines(x_seq,vel1bis,lwd=2,col="red",type="b")
+#lines(x_seq,vel2,lwd=2,col="black")
+lines(x_seq,vel2bis,lwd=2,col="black",type="b")
+legend(0.4, 1, c(expression("Numerics"), expression("Alipchenkov eq 6"), expression("Alipchenkov eq 15")), 
+lwd= c(2, 2,2),col=c("blue","red","black" ))
+
+dev.off()
+######################################################################
+######################################################################
+
+pdf("comparison-flux-wrong.pdf")
+plot(x_seq,flux_v,col="blue",xlab=expression("x"),ylab=expression(" average flux"),
+ylim=range(c(min(flux_v,flux),max(flux_v,flux))))
+lines(x_seq,flux_v,col="blue", type="b")
+lines(x_seq,flux,col="red")
+lines(x_seq,flux,col="red",type="b")
+dev.off()
+
+
+pdf("comparison-flux-correct.pdf")
+plot(x_seq,flux_v_correct,col="blue",xlab=expression("x"),ylab=expression(" average flux"),
+ylim=range(c(min(flux_v_correct,flux_correct),max(flux_v_correct,flux_correct))),type="b")
+#lines(x_seq,flux_v_correct,col="blue", type="b")
+#lines(x_seq,flux_correct,col="red")
+lines(x_seq,flux_correct,col="red",type="b")
+dev.off()
+
+
+
+pdf("comparison-density.pdf")
+plot(x_seq,dens_v,col="blue",xlab=expression("x"),ylab=expression(" density"),
+ylim=range(c(min(dens_v,dens),max(dens_v,dens))),type="b")
+lines(x_seq,ana_dens,lwd=2,col="black")
+lines(x_seq,dens,col="red",type="b")
+dev.off()
+
+
+
+pdf("comparison-sig.pdf")
+plot(x_seq,sig_v,col="blue",xlab=expression("x"),ylab=expression(" average sig_vv squared"),
+ylim=range(c(min(sig_v,sig),max(sig_v,sig))),type="b")
+#lines(x_seq,sig_v,col="blue", type="b")
+#lines(x_seq,sig,col="red")
+lines(x_seq,sig,col="red",type="b")
+dev.off()
+
+pdf("comparison-vel-squared.pdf")
+plot(x_seq,velsq_v,col="blue",xlab=expression("x"),ylab=expression(" average v squared"),
+ylim=range(c(min(velsq_v,velsq),max(velsq_v,velsq))),type="b")
+#lines(x_seq,velsq_v,col="blue", type="b")
+#lines(x_seq,velsq,col="red")
+lines(x_seq,velsq,col="red",type="b")
+dev.off()
+
+pdf("comparison-vel.pdf")
+plot(x_seq,vel_v,col="blue",xlab=expression("x"),ylab=expression(" average v"),
+ylim=range(c(min(vel_v,vel),max(vel_v,vel))),type="b")
+#lines(x_seq,vel_v,col="blue", type="b")
+#lines(x_seq,vel,col="red")
+lines(x_seq,vel,col="red",type="b")
+dev.off()
+
+
+pdf("comparison-vel-bis.pdf")
+plot(x_seq,vel_v,col="blue",xlab=expression("x"),ylab=expression(" average v"),
+ylim=range(c(min(vel_v,vel),max(vel_v,vel))),type="b")
+lines(x_seq,0.2*x_seq,col="green")
+#lines(x_seq,vel,col="red")
+lines(x_seq,vel,col="red",type="b")
+dev.off()
+
+myseq2<-mymat_save[ ,2]
+norm1<-def_integral(myseq2,h)
+
+myseq2<-mymat[ ,2]
+norm2<-def_integral(myseq2,h)
+
+pdf("comparison-vel-profiles.pdf")
+plot(mymat_save[ ,1],mymat_save[ ,2]/norm1
+,col="blue","l",xlab=expression("v"),ylab=expression("f(x,v,t) at wall"))
+lines(mymat[ ,1],mymat[ ,2]/norm2,col="red")
+dev.off()
+
+
+
+
+print("So far so good")
\ No newline at end of file