スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

12月2日

12月2日、ヒストグラム作成時に使用した全ファイル:

test3.c

#include (←studio.hがかぎ括弧付きで書いてある)
#include (←math.hがかぎ括弧付きで書いてある)

#define RPNUM 20
#define MAX 400000

int main(void){
int i,j,k,target;
FILE *fpi;
FILE *fpG;
FILE *fpV;
FILE *fpT;
char filenamei[50]; //書込み用
char filenameG[50];
char filenameV[50];
char filenameT[50];

double gyr,theta,V;
int gen;

for(i=0; i<1; i++){
for(target=12; target<=12; target++){
gen=0;
sprintf(filenameG,"CNMD_G_%d.dat",i);
sprintf(filenameV,"CNMD_V_%d.dat",i);
sprintf(filenameT,"CNMD_T_%d.dat",i);
sprintf(filenamei,"hist2dCNMD%d_T%d.dat",i,target);

fpG=fopen(filenameG,"r");
fpV=fopen(filenameV,"r");
fpT=fopen(filenameT,"r");
fpi=fopen(filenamei,"w");

for(k=0; k fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
}
printf("%lf %lf %lf\n",gyr,V,theta); //TARGETで指定したRPの温度を確認する

for(k=0; k<(int)RPNUM-target; k++){ //温度の項を脱出
fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
}

for(k=0; k fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
}

for(j=0; j<100000; j++){
for(k=0; k<(int)RPNUM-1; k++){
fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
}
fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
gen +=10;
}

for(j=0; j for(k=0; k<(int)RPNUM-1; k++){
fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
}
fscanf(fpG,"%lf",&gyr);
fscanf(fpV,"%lf",&V);
fscanf(fpT,"%lf",&theta);
if((gyr>200 && gyr<550)&&(theta>20 && theta<30))
fprintf(fpi,"%d %lf %lf %lf\n",gen,gyr*10.0,V,theta*100.0);
gen +=10;
}

fprintf(fpi, "\n");
for(k=0; k<32000; k++){ //θ固定、gyrをプロット
gyr = 2000.0 + (double)k/32000.0*3500.0;
fprintf(fpi, "%d %lf 0 2899.0\n", -k, gyr);
}

for(k=0; k<40000; k++){ //θ固定、gyrをプロット
theta = 2000.0 + (double)k/40000.0*1000.0;
fprintf(fpi, "%d 4999.0 0 %lf\n", -32000 - k, theta);
}

fclose(fpG);
fclose(fpV);
fclose(fpT);
fclose(fpi);
}
}
return 0;
}



このプログラムを実行して、元々の巨大データを割って独立した時系列データに換え、Rにわたす。範囲を指定するマイナス値の時間の部分(追加で書き足す時系列部分)で指定する領域は、Rのhist2d()にわたすことを考えると、binsが計算時間上固定されているために格子が荒くなるという意味で狭い方が本当はいいのだが、時間もないので広くとることにした。

install.Rは前の記事で書いてある。

test3.R

data1 <- read.table("hist2dCNMD0_T12.dat")
datax <- data1[,2]
datay <- data1[,4]
h2d <- hist2d(datax, datay, show=FALSE, same.scale=TRUE, nbins = c(800,800))
zzz <- h2d$counts
xxx <- h2d$x
yyy <- h2d$y
dataX <- matrix(nrow=640800,ncol=3)
for(i in 1:800){
for(j in 1:800){
dataX[i+(800*j),1] <- xxx[i]
dataX[i+(800*j),2] <- yyy[j]
if(zzz[i,j] == 0)
dataX[i+(800*j),3] <- 1
else{
if(zzz[i,j] >= 250)
dataX[i+(800*j),3] <- 250
else
dataX[i+(800*j),3] <- zzz[i,j]
}
}
}
data3 <- as.data.frame(dataX)



その後は手動で以下を実行する(12/2に作成したヒストグラムは以下のコマンドで作成した):

sun <- c("#110000","#660011","#bb0000","#C93211","#D76423","#F0BB43","#ffee55","#FFF5A3","#FFFCE8")
sea <- c("#002046","#0066ff","#11aadd","#44ccee","#00ddff","#77ffff","#bbeeee","#ffffee")
alien <- c("#000000","#447744","#559944","#55dd33","#bbcc55","#ffff77","#ffffff")
darkmoon <- c("#000011","#336699","#5599bb","#aaccdd","#EBEDD4","#F3F3E7","#eeffff")


ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, fill=data3[,3]/400000.0)) + geom_tile() + xlim(c(20,32)) + ylim(c(21.0,28.9))+ opts(title="T=0.60") + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + scale_fill_gradientn("frequency",limits = c(0, 0.000625) , colours = sea)



darkmoon

Rでの様々な記号の出力法を知った。今回は無限記号をタイトルに出力した。

> ggplot(data3, aes(x=data3[,1]/150.0, y=data3[,2]/100.0, fill=data3[,3]/400000.0)) + geom_tile() + xlim(c(20,32)) + ylim(c(21.0,28.9))+ opts(title=expression(paste("CageSize=",infinity,", T=0.65"))) + xlab(expression(R[g]^2)) + ylab(expression(paste(Sigma,theta))) + scale_fill_gradientn("frequency",limits = c(0, 0.000625) , colours = sea)


NMD_65_2_convert_20101203021620.png


参考にしたサイト
mathematical annotation in R

また、今日参考にしたサイト

texの数式中の括弧の大きさの調整:\right[ \left] とすると自動で大きさが変わる
Bracket: TeX

volcanoのデータが綺麗に3Dで描かれている:
R-Source RGL
スポンサーサイト

コメントの投稿


非公開コメント

コメント

Profile

zoa

Author:zoa
Metropolis (1927)

Calendar
05 | 2017/06 | 07
- - - - 1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 -
Labels
Admin
Previous Posts
Recent Comments
Recent Trackback
Archive
このページのトップへ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。