制限三体問題シミュレータ by Rust

// 制限三体問題シミュレータ
//
// 座標系は連星系の重心を原点とし, 平均運動nで回転する系.
// 従ってe=0のとき連星は静止するが, e \neq 0も計算できる.
// 連星の長半径aと平均運動nが1となる単位系を用いる.
//
// 必要なパラメータは
//     連星の質量比 q = M_2 / M_1
//     連星の離心率 e
//     計算終了時刻 t_end
// と初期条件 x: [f64;3], v: [f64;3] の計九つ. 上ふたつを
// コマンドライン引数としてセットすることもできる.
//

// use std::env; // コマンドライン引数を使用する場合はコメントアウトを外す

use std::f64::consts::PI;

fn main() {
 // パラメータの設定
 // let args: Vec = env::args().collect();
 // let q: f64 = args[1].parse::().unwrap();
 // let e: f64 = args[1].parse::().unwrap();
 let q: f64 = 0.01;
 let e: f64 = 0.;
 let t_end: f64 = 128.*PI;

 // 初期条件
 let mut t: f64 = 0.;
 let mut x: [f64;3] = [ 0.5-q/(1.+q), 3f64.sqrt()/2., 0. ];
 let mut v: [f64;3] = [ 0., 0., 0. ];
 println!("{} {} {} {}", t, x[0], x[1], x[2] );

 // アウトプット間隔
 let t_output: f64 = 0.0625;
 let mut t_record: f64 = t_output;

 // 数値積分
 while t < t_end {
  rk4( &mut t, &mut x, &mut v, [q,e] );
  if t >= t_record {
   println!("{} {} {} {}", t, x[0], x[1], x[2] );
   t_record += t_output;
  };
 };
}


fn true_anomaly ( t: f64, e: f64 ) -> f64 {
 let mut ecc: f64 = t;
 let mut delta_e: f64 = - ( ecc - e*ecc.sin() - t ) / ( 1. - e*ecc.cos() );
 while delta_e.abs() > 1e-4 {
  ecc += delta_e;
  delta_e = - ( ecc - e*ecc.sin() - t ) / ( 1. - e*ecc.cos() );
 };
 ( ((1.+e)/(1.-e)).sqrt() * (ecc/2.).tan() ).atan() * 2.
}


fn force( t: f64, x: [f64;3], v: [f64;3], param: [f64;2] ) -> [f64;3] {
 // param: [f64;2] = [ q, e ];
 let x1: [f64;3];
 let x2: [f64;3];
 { // 連星の位置を計算
  let f: f64 = true_anomaly( t, param[1] );
  let r: f64 = (1.-param[1].powi(2))/(1.+param[1]*f.cos());
  x1 = [ -param[0]/(1.+param[0])*r*(f-t).cos(), -param[0]/(1.+param[0])*r*(f-t).sin(), 0. ];
  x2 = [ 1./(1.+param[0])*r*(f-t).cos(), 1./(1.+param[0])*r*(f-t).sin(), 0. ];
 }
 let mut f: [f64;3] = [ x[0]+2.*v[1], x[1]-2.*v[0], 0. ];
 { // x1とx2からの力
  let r1: f64 = 1./(1.+param[0])*( (x[0]-x1[0]).powi(2) + (x[1]-x1[1]).powi(2) + (x[2]-x1[2]).powi(2) ).powf(-1.5);
  for k in 0..3 { f[k] += - r1*(x[k]-x1[k]); }
  let r2: f64 = param[0]/(1.+param[0])*( (x[0]-x2[0]).powi(2) + (x[1]-x2[1]).powi(2) + (x[2]-x2[2]).powi(2) ).powf(-1.5);
  for k in 0..3 { f[k] += - r2*(x[k]-x2[k]); }
 }
 f
}


fn rk4( t: &mut f64, x: &mut [f64;3], v: &mut [f64;3], param: [f64;2] ) {
 let x1: [f64;3] = *v;
 let v1: [f64;3] = force( *t, *x, *v, param );

 let h: f64 = 0.0625*((1f64.min(1./v1[0].abs())).min(1./v1[1].abs())).min(1./v1[2].abs());

 let x2: [f64;3] = [ v[0]+v1[0]*h*0.5, v[1]+v1[1]*h*0.5, v[2]+v1[2]*h*0.5 ];
 let v2: [f64;3] = {
  let xp: [f64;3] = [ x[0]+x1[0]*h*0.5, x[1]+x1[1]*h*0.5, x[2]+x1[2]*h*0.5 ];
  force( *t+h*0.5, xp, x2, param )
 };

 let x3: [f64;3] = [ v[0]+v2[0]*h*0.5, v[1]+v2[1]*h*0.5, v[2]+v2[2]*h*0.5 ];
 let v3: [f64;3] = {
  let xp: [f64;3] = [ x[0]+x2[0]*h*0.5, x[1]+x2[1]*h*0.5, x[2]+x2[2]*h*0.5 ];
  force( *t+h*0.5, xp, x3, param )
 };

 let x4: [f64;3] = [ v[0]+v3[0]*h, v[1]+v3[1]*h, v[2]+v3[2]*h ];
 let v4: [f64;3] = {
  let xp: [f64;3] = [ x[0]+x3[0]*h, x[1]+x3[1]*h, x[2]+x3[2]*h ];
  force( *t+h*0.5, xp, x4, param )
 };

 *t += h;
 for k in 0..3 {
  x[k] += ( x1[k] + 2.*( x2[k] + x3[k] ) + x4[k] ) * h / 6.;
  v[k] += ( v1[k] + 2.*( v2[k] + v3[k] ) + v4[k] ) * h / 6.;
 }
}