open その２

テストプログラムを作って見ました

これは，振り子が動く図をたくさん出力します
これをimagemagickのconvertを使って動画にします
たとえば，

convert -delay 10 IMGSEQUENCE/*.gif a.mpg

こんな感じです

ちょっと複雑になりました

#include ＜stdio.h＞
#include ＜stdlib.h＞
#include ＜math.h＞
#include ＜string.h＞

int main(void) {
  FILE *gid, *data;
char *data_file, *k,kk[20],k1[10],k2[10],k3[10];
float dx, x, y,th,vth,al,dt;
int i, imax,j ,jmax,j1,j2,j3;
 imax=1;
 jmax=200;
  gid = popen("gnuplot -persist", "w");

    fprintf(gid, "set xrange [-2:2]\n");
fprintf(gid, "set yrange [-2:2]\n");
    fprintf(gid, "set term gif truecolor size 640,480\n");
data_file="out.dat";
  
  
     dx=0.4/imax;
dt=0.2;
al=0.1;
    x=0;
th=1;
vth=0;
for(j=0; j
j1=j%10;
j2=((j-j1)%100)/10;
j3=j/100;
sprintf(k1,"%d",j1);
sprintf(k2,"%d",j2);
sprintf(k3,"%d",j3);


sprintf(kk,"data/outdata%d",j);

data_file=kk;
    data = fopen(data_file,"w");
fprintf(data,"%f\t%f\n", 0.0 ,0.0);
 //    for(i=0; i<=imax; i++){
      th+=vth*dt;
 vth+=-al*th*dt;
 x=sin(th);
 y=-cos(th);
       fprintf(data,"%f\t%f\n",x, y);
       x+=dx;
//     }
 
     fclose(data);
   sprintf(k,"%d",i);
 //   printf("%s\n",k);
strcat(k,k3);
strcat(k,k2);
strcat(k,k1);
    fprintf(gid, "set output 'IMGSEQUENCE/sequence0%s.gif' \n",k);

    fprintf(gid,  " plot '%s' w l \n",kk);
    fprintf(gid, "set output \n");
};
  pclose(gid);
  return 0;
}

popen

popen・・・Gnuplotを用いてアニメを作成するときに使うC言語の関数
各自調べて勉強してみましょう！！

ポインタ・数値微分についてと次回への課題

ポインタ
ポインタを使ってファイルに書き込む作業（過去の記事参照）
ポインタ：メモリのアドレス
イメージはファイルの置いてある場所
ex　棚の右から3番目とか

数値微分
微分を差分に置き換える。
微分はf(x,t)=dx/dt
差分はxi+1-xi=f(x,t)Δt
     詳しくは過去記事に数式が記載されている。

次回の課題
二次元振り子（近似なし）を差分プログラムを作成。

昨日のプログラムは動きました

以下のプログラムは、ぼくの計算機ではコンパイルと実行ができました
ので、昨日のプログラムには、文字化けがあるような気がします
#include のあとが問題かもしれません
羽部

ファイルへの書き込み

＜今日の内容＞
・sinxのテーラー展開の復習
・得られたデータをファイルに書き込む
・Gnuplotでデータをプロットする

〜ファイルへの書き込み〜
「./a.out > ”ファイル名”」とすればリダイレクトされて、./a.outの内容がファイルに書き込まれる。

sinxの展開のプログラムでpow関数を用いた方がコンパイルできなかった。
→各自、原因を考えて来てください

プログラム
#include
#include
main()
{
int i, n;
double x;
double dx, s ,c1,c3,c5,c7,c9;
x=0.0e0;
n= 70;
dx = 0.1e0;
x = -dx ;
c1=1;
c3=c1/2/3;
c5=c3 / 4 /5;
c7 = c5 /6 /7;
c9 = c7 /8 /9;
for (i=1;i<=n; i++ ){
x = x + dx;
s = x -c3 *pow(x,3.0) + c5*pow(x,5.0) - c7 * pow(x,7.0)+c9* pow(x,9.0);
printf(" %d %f %e \n", i , x, s);
}
}

前回の続き

・ドリフト運動
　ドリフト運動について電荷にかかる合力や曲率半径などを考えて、ドリフト運動の理解を深めようとした。
　曲率半径について
　　ds = 2πR *（dθ／２π）= R＊dθ
　
　　R * |de/ds| = 1 (e:速度の単位ベクトル)←相似比より

　　|de/ds| * ds = dθ←上の２式でRを消去


・ドリフト運動（その２）
　x方向の初速を０にした場合をグラフで表した。

・ドリフト運動（その３）
　x方向の初速を−２．０にした場合をグラフで表した。

ドリフト運動を調べる（概要）

まず電場がかかっていない、普通のローレンツ力の場合を考えて比較する。

ローレンツ力は円運動
なぜなら速度と向心力の内積が０←向心力は速度に影響しない

円運動とローレンツ力のつりあいの式
mv^2/r = qvB
より
r=m/(qB)×v
等速円運動（r,vは定数）

ドリフト運動は向心力の他に電場による静電気力も働くので内積は０にならない。（←あとで証明）

速度vは時刻とともに変化する。←変化は電場による加速度の速度方向をみて

r=m/(qB)×v

より、速度が大きくなれば荷電粒子の回転半径が大きくなることが分かる。

gnuplotで各軸の感覚を調節する

グラフが正しい比率でかかれないと
物体の運動の軌跡が解析しにくいことがある。

これにはグラフで表示するxとyの範囲を適度に調節すればよい。
x軸の範囲を0以上0.6以下にしたいときは、
以下のコマンドをプロット前に実行しておく。

set xrange [ 0 : 0.6 ]

y軸についても同様に行えばよい。

gnuplotで各軸に名前をつける

gnuplotは標準では各軸にX軸、Y軸などの名前が表示されない。
gnuplotでデータをプロットする前に、次のコマンドを実行しておく。

set xlabel "x" 0.0,1.0

これでx軸の下に”x”と表示されるようになる。

ドリフト運動のグラフを1:1に直してみると…

磁場の加速度のx方向の成分が電場の加速度の方向と強めあうときは、グラフの接線方向の加速度が大きくなるので回転半径が大きくなる。また弱め合う場合は、グラフの接線方向の加速度が小さくなるので回転半径が小さくなる。この変化を繰り替えしていくと荷電粒子が全体としてy軸負の方向に移動しているようにみえる。

各時刻でグラフに内接する円の中心の座標がどのように移動していくかプロットすることで、何か分かるだろうか…？

ドリフト運動（電場×磁場;垂直）（重力あり）

本日のソースコード

#include <>

main(){

double a,B,E,g,m,q,C,x,y,z,Vx,Vy,Vz,Vx1,Vy1,dt;

//C=q*B/m（定数）

int i,n;

a=0.0; //空気抵抗
B=2.0; //磁束密度//z方向
E=1.0; //電場//x方向
g=9.8; //重力加速度
m=1.0; //粒子質量
q=5.0; //粒子電荷
x=0.0; //
y=0.0; //
z=0.0; //
Vx=2.0; //
Vy=0.0; //
Vz=0.0; //


dt=0.0001; //時間（パラメータ）
for(i=0;i<100000;i++)
{
C=q*B/m;
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt+q*E*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
Vz=Vz-g*dt-a*Vz*dt;
printf("%e %e %e \n", x, y, z );
}
}

x-y平面への出力結果
x方向へ少しずつ広がっているのは微小時間tの値が大きいことによる誤差と考えられます。

x軸+方向に電場をかけた結果、ドリフト運動してy軸-方向へドリフト運動します。
さらに、振動の中心がすこしx軸+方向へずれているのが分かります。

z軸-方向へ重力場がはたらいているのでだんだん上昇速度が小さくなっています。

三次元plotその2 空気抵抗あり

右の図は下のプログラムの場合です


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

三次元のplot

先週のグラフを三次元で表してみました
（z軸方向に等速度運動）

gnuplot > plot "re"
→gnuplot > splot "re"
で三次元のplotが表示される
splotはファイルreに保存されたデータから一行目をx座標、二行目をy座標、三行目をz座標としてplotする。


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.0;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

来週までの宿題

・プロットする点の減らし方
・磁場が変化する場合
・"三"次元プロット

磁場中の荷電粒子の運動（空気抵抗あり）

一様磁場：画面手前方向
電荷：＋
初速度：＋x軸方向

#include

main(){

double a,x,y,Vx,Vy,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
Vx=2.0;
Vy=0.0;
C=10;
dt=0.0001;
for(i=0;i<1000000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;

Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e \n", x, y);
}

2009年10月26日月曜日へのコメント

return 0は関数の実行結果の値をゼロにするという意味だと思います
main関数でこうしてもあまり意味がないように思います

stdio.hは標準の入出力で利用する関数を呼び出すために必要な
情報を与えていると思えば良いと思います

来週の予告

磁場中の荷電粒子の運動（ローレンツ力あり）

運動方程式（xとyの連立方程式）
m*(dVx/dt)= q*Vy*B
m*(dVy/dt)=-q*Vx*B

どうやって解くか考えましょう。

重力場中（空気抵抗あり）の物体のx-tグラフのための数値計算

#include < stdio.h >
main()
{
double t, x, dt, dx, ta, tb, v;
double f(double ,double );
int i, n;
x = 0.0 ;
t = 0.0 ;
ta = 0.0 ;
tb = 100 ;
v = 0.0 ;
n = 1000;
dt = (tb - ta) / n ;
for ( i = 0; i < n ; i++){
x = x+ v*dt ;
v = v+ f(v, t)*dt ;
t = t + dt ;
printf(" %e %e \n", t, x);
}
}
double f( double x, double t){
double y, z, alpha, g;
alpha = 0.1;
g = - 9.8;
z = -alpha * x + g ;
return z;
}

来週のメニュー

来週は
2階の微分方程式（1次の項を含む）を
電解中の荷電粒子の運動（空気抵抗あり）をモデルに
やります。

ヒント
for()
{
x = x + v*dt
v = v + (A-B*v)*dt

}

(A=qE/m, B=alpha/m)

グラフを画像ファイルとして保存する

gnuplotで表示した画像の保存
clipboardに保存（グラフ画面の左上のアイコンをクリック）
GIMPを起動しファイルメニューから”画像の生成”を選び”クリップボードから”をクリック。するとグラフが表示される。名前をつけて保存を選び、形式を”JEPG画像”にして保存。

ブログにアップロードするには、ブログの投稿画面で”画像を追加”から”参照”