/* ファイル名はkinema.c */
#include <stdio.h>
#include <math.h>

double calcE(double M, double p);
double calcP(double M, double E);

kinema()  /*  <--- main文の名前をかえる   */
{
  double M1, M2, M3, M4;  /* 粒子の質量 */
  double p1, p2, p3, p4;  /* 粒子の実験室系での運動量 */     
  double E1, E2, E3, E4;  /* 粒子の実験室系でのエネルギー */ 
  double theta;           /* 実験室系での散乱角 */           
  double p3_parl;         /* 粒子３の実験室系での運動量の
			     重心系の速度に射影した成分*/

  double p_cm1, p_cm2, p_cm3, p_cm4; /* 粒子の重心系での運動量 */     
  double E_cm1, E_cm2, E_cm3, E_cm4; /* 粒子の重心系でのエネルギー */ 
  double theta_cm;		     /* 重心系での散乱角 */           

  double beta_cm, gamma_cm;          /* 重心系の速度 */
  double E_cm;
			     
  /* input */
  M1 =   938.2;
  M2 = 11177.0;
  M3 =   938.2;
  M4 = 11177.0;


  p1 = calcP(M1, M1+12.0);   /* 運動エネルギー12MeV */
  p2 = 0.0;

  E1 = calcE(M1, p1);  /* 質量と運動量からエネルギーを求める関数 */
  E2 = calcE(M2, p2);   

  
  /* output fileを開く */
  FILE *fp;
  fp = fopen("output.txt","w");
    
  /* まずbeta, gammaを求める */
  beta_cm = (p1+p2)/(E1+E2);
  gamma_cm = 1.0/sqrt(1.0-beta_cm*beta_cm);

  /* 重心系のエネルギーを求める */
  E_cm = sqrt((E1+E2)*(E1+E2)-(p1+p2)*(p1+p2));

  /* 散乱粒子の重心系でのエネルギーを求める */
  E_cm3 = (E_cm + (M3*M3 - M4*M4)/E_cm)/2.0;
  E_cm4 = (E_cm - (M3*M3 - M4*M4)/E_cm)/2.0;
  p_cm3 = p_cm4 = calcP(M3, E_cm3); /* 質量とエネルギーから運動量を求める関数 */
  for (theta_cm=0.0; theta_cm<=180.0; theta_cm+=1.0) {
    /* 重心系でのエネルギー、運動量に戻そう */
    p3 = sqrt(pow(gamma_cm*beta_cm*E_cm3+gamma_cm*p_cm3*cos(theta_cm*3.14/180.0),2.0) + 
	      pow(p_cm3*sin(theta_cm*3.14/180.0), 2.0));

    p3_parl = gamma_cm*beta_cm*E_cm3+gamma_cm*p_cm3*cos(theta_cm*3.14/180.0);
    if (p3_parl >= 0)
      theta = asin((p_cm3/p3)*sin(theta_cm*3.14/180.0))*180.0/3.14;
    else
      theta = 180.0 - asin((p_cm3/p3)*sin(theta_cm*3.14/180.0))*180.0/3.14;

    printf("%f %f\n", theta, calcE(M3, p3)-M3); /*  <--- 実験室系の角度と運動エネルギーを表示  */

  
    /* output fileに書き出す */	 
    fprintf(fp, "%f %f\n", theta, calcE(M3, p3)-M3);
  
   }
  
    /* output fileをとじる */	 
    fclose(fp);
  
}

double calcE(double M, double p)
{
  return sqrt(M*M + p*p);
}

double calcP(double M, double E)
{
  return sqrt(E*E -M*M);
}