###############################################################
# beams.awk                                                   #
# 球形の風船によって屈折させられる波の進路のシミュレーション  #
# (C)1997.  Hideo Konami, MNCT                                #
###############################################################
#
#              計算の手順
#
# 1. 計算のための次のパラメータを適宜設定する
#    物理的パラメータ:  R, M1, M2
#    描画用: rR, wx, wy, fx, fy, nbeam
#
# 2. ビームの形状を決める関数を選択する。
#
# 3. プログラムを実行して、LaTeX の picture 環境のフォーマット
#    のファイルを生成する。
#    awk -f beams.awk > hoge.tpc
#
# 4. eepic.sty を取り込んだ LaTeX ファイル中に hoge.tpc を input
#    しておく。     (\input hoge.tpc)
# 
# 5. LaTeX ファイルをコンパイルして、結果を表示する。
#

BEGIN{
  Pi = 3.14159265358979;
  Pi180 = 180/Pi;

# パラメータ設定
# 風船の物理半径
  R = .5;
# M1: 外側の気体の分子量, M2: 内側の気体の分子量  
  M1 = 29.0;
  M2 = 44.0;
# 風船の直径と図の縦サイズの比
  rR = 0.6;
# wx: 横(mm), wy: 縦(mm) 描画する枠のサイズ
  wx = 170;
  wy = 100;
# 球の中心を (wx*fx,wy*fy) の位置に置く
  fx = 0.6;
  fy = 0.5;
# ビームの本数
  nbeam=51;

  SetParameters();

# ビームを生成する。
# 下の２つの定義関数のどちらかを使って走らせること。
# '#' が先頭にあると、その行は実行されないので、実行させたい
#  方の行の先頭の '#' を取り除き、そうでない方は付けておく。

#　平行光束を生成する。 
#  GeneratePararellBeams(i,j,k);

# 一点から広がる光束を生成する。このときは
# 光源の位置を示す distance も設定しなければならない。
# distance: 球の中心から半径の何倍のところにあるか。
  distance = R * 1.95;
  GenerateSpreadBeams(i,distance);

#    　　　　　　設定するのはここまで
#-----------------------------------------------------------

# 球面を描く
  DrawSphere();
  
# ビームの経路の計算
  for ( i=1;i<=nbeam;i++){

#   判別式を調べて、正なら計算を行う。
    if ((D = Det(i)) > 0) {
#     (x1,y1): 入射する点の座標
      x1 = (-a1[i]*b1[i]+sqrt(D))/(a1[i]^2+1);
      y1 = a1[i]*x1+b1[i];

#     th1: 入射角 
      th1 = atan2(a1[i],1)-atan2(y1,x1);
#     eta1: 出射角,  t: 軸と出射光線のなす角
      se = sin(th1)/rho;
      if (se^2 > 1 ){
	DrawBeam1(i);
	continue;
      }
      ce = sqrt(1-se^2);
      eta1 = atan2(se,ce);

      t = atan2(y1,x1) + eta1;
      a2 = sin(t)/cos(t);
      b2 = y1 - a2 * x1;

      x2 = -(x1+a2*y1+a2*b2)/(1+a2^2);
      y2 = a2 * x2 + b2;


#   風船から出て行くビーム
#   球面 C 上のP2:(x2,y2) に傾き a2 で入射するビームの出ていく
#   方向を調べる。
      th2 = atan2(a2,1) - atan2(y2,x2); 
      se = sin(th2) * rho;
      if(se^2 >1){
	DrawBeam2();
	continue;
      }

      ce = sqrt(1 - se^2);
      eta2 = atan2(se,ce);
      t = atan2(y2,x2) - eta2;
      a3 = sin(t)/cos(t);
      b3 = y2 - a3 * x2;
#      print "theta2,eta2 = ", th2, eta2;

#   線を描画する。
      DrawBeam1(i);
      DrawBeam2();
      DrawBeam3();
    }
    else{
      DrawSingleBeam(i);
    }
  }
  printf("\\end{picture}\n");
}
function DrawSingleBeam(i){
  printf("%% Non-crossing beam a1[%d] = %.2f\n",i,a1[i] );
  a = a1[i];
  b = b1[i];
  wxr = wx * ax + bx;
  wxl = bx;
  wyu = wy * ay + by;
  wyl = by;
  #print ">>>>"  b * Ay + By;
  if (a >= 0) {
    if ((yr = a * wxr + b) <= wyu){
      xr = wxr;
    }
    else{
      yr = wyu;
      xr = (wyu-b)/a;
    }
    if((yl = a * wxl + b) > wyl){
      xl = wxl;
    }
    else{

      yl = wyl;
      xl = (wyl-b)/a;
    }
  }
  else{
    if((yr = a * wxr + b) >= 0){
      xr = wxr;
    }
    else{
      yr = 0;
      xr = -b/a;
    }
    if((yl = a * wyl + b) <= wyu){
      xl = 0;
    }
    else{
      yl = wyu;
      xl = (wyu-b)/a;
    }
  }
  printf("\\put(0,0){\\path(%.2f,%.2f)(%.2f,%.2f)}\n",
	 xr*Ax+Bx,yr*Ay+By,
	 xl*Ax+Bx,yl*Ay+By);
}

function DrawBeam1(i){
  a = a1[i];
  b = b1[i];
  wxr = wx * ax + bx;
  wyu = wy * ay + by;
  wyl = 0 * ay+ by;
  printf("%% DrawBeam1\n");
  printf("\\allinethickness{0.4mm}\n");

  if (a >= 0){
    if ((yt = a *wxr + b)<=wyu){
      xt = wxr;
    }
    else{
      yt = wyu;
      xt = (wyu-b)/a;
    }
  }
  else{
    if((yt = a * wxr + b) >= wyl){
      xt = wxr;
    }
    else{
      yt = wyl;
      xt = (wyl-b)/a;
    }
  }
  printf("\\put(0,0){\\path(%.2f,%.2f)(%.2f,%.2f)}\n",
	 x1*Ax+Bx,y1*Ay+By,
	 xt*Ax+Bx,yt*Ay+By);
}
function DrawBeam2(){
  printf("%% DrawBeam2\n");
  printf("\\put(0,0){\\path(%.2f,%.2f)(%.2f,%.2f)}\n",
	 x1*Ax+Bx,y1*Ay+By,
	 x2*Ax+Bx,y2*Ay+By);
}
function DrawBeam3(){
  a = a3;
  b = b3;
  wxl = bx;
  wyu = ay * wy + by;
  wyl = by;
  
  yl0 = a3 * wxl + b3;
  if(yl0 >= wyl && yl0 <= wyu){
    xl = bx;
    yl = yl0;
  }
  else if (yl0 > wyu){
    xl = (wyu-b3)/a3;
    yl = wyu;
  }
  else{
    xl = (wyl-b3)/a3;
    yl = wyl;
  }
  printf("%% DrawBeam3\n");
  printf("\\put(0,0){\\path(%.2f,%.2f)(%.2f,%.2f)}\n",
	 x2*Ax+Bx,y2*Ay+By,
	 xl*Ax+Bx,yl*Ay+By);
}


function GeneratePararellBeams(i,j,k){
  method = "GeneratePararellBeams";
  a0 = 0.0;
  db = wy * ay /(nbeam+1);
  bs = db + by;

  for (i=1;i<=nbeam;i++){
    a1[i] = a0;
    b1[i] = (i-1)*db + bs;
    #print b1[i];
    #print ">>" b1[i]*Ay + By;
  }
}
function GenerateSpreadBeams(i,distance){
  method = "GenerateSpreadBeams";
  wyu = ay * wy + by;
  th  = 2*atan2(wyu,distance);
  dth = th/(nbeam-1);
  th0 = -th/2;

  for (i=1;i<=nbeam;i++){
    th = th0 + (i-1)*dth;
    a1[i] = sin(th)/cos(th);
    b1[i] = -a1[i]*distance;
    #print b1[i];
    #print ">>" b1[i]*Ay + By;
  }
}

function Det(i){
  return  a1[i]^2 * b1[i]^2 - (a1[i]^2 + 1)*(b1[i]^2-R^2);
}

# ヘッダーと円、いくつかのコメントを出力する。
function DrawSphere(){
  pline = "%%%%%%%%%%%%%%%%%%%%%%%%";
  printf("%s%s\n",pline,pline);
  printf("%%       Beam  Refraction by a Spheric Lens     %%\n");  
  printf("%% Generated with beams.awk.(C)1997 Hideo Konami%%\n");  
  printf("%s%s\n%%\n",pline,pline);

  printf("%%  == Parameters in this calculation ==  \n");
  printf("%%     rR = %.2 , nbeam = %d\n",rR,nbeam);
  printf("%%     wx = %.2 mm, wy = %.2 mm\n",wx,wy);
  printf("%%     fx = %.2 mm, fy = %.2 mm\n",fx,fy);
  printf("%%     method = '%s'\n",method);
  if (method = "GenerateSpreadBeams"){
    printf("%%     distance = %.2f\n",distance);
  }
  printf("\\unitlength1mm\n");
  printf("\\begin{picture}(%.2f,%.2f)(0,0)\n",wx,wy);
  printf("\\put(0,%.2f){\\makebox(%.2f,5)",wy+1,wx);
  printf("{$M_1 = %.2f, M_2 = %.2f\\ (v_1/v_2 = %.4f)$}}",
	 M1,M2,rho);
  printf("\\put(0,-5){\\makebox(%.2f,4)[r]{\\footnotesize\\bf",wx);
  printf(\
  " wx = %.2f, wy = %.2f, Radius = %.2f, center =(%.2f,%.2f) }}\n",
	 wx,wy,rR*wy/2,wx*fx,wy*fy);

  printf("\\allinethickness{0.6mm}\n");
  printf("\\put(0,0){\\framebox(%.2f,%.2f){}}\n",wx,wy);

  printf("\\allinethickness{0.3mm}\n");
  printf("\\put(%.2f,%.2f){\\circle{%.2f}}\n",cx,cy,R*Ax*2);

}

function SetParameters(){

  rho = sqrt(M2/M1);
  cx = wx * fx;
  cy = wy * fy;
  wR = wy * rR /2;

# Drawing scale to physical space.
  Ax = wR/R;
  Bx = cx;
  Ay = wR/R;
  By = cy;
# Physical scale to drawing scale.
  ax = R/wR;
  bx = -ax*cx;
  ay = R/wR;
  by = -ay*cy;
}