#define q 1.60217653e-19  /* coulomb */
#define me 0.510998902 /* e MeV/c^2 */
#define mp 938.271998  /* ion MeV/c^2 */
#define ma 3727.409 /* alpha MeV/c^2 */
#define c 2.99792458e8 /* m/s */
#define B0 1. /* arbitary unit */
#define AU 1.496e11 /* 1 AU in meter */
#define lperp 0.02 /* bend over length scale (AU)*/
#define lz 0.00167 /* AU */ 
#define B00 5.e-9  /*tesla -- for nondimensionalization */
#define pi 4.*atan(1.0) 

main()
{
   double E,KE,m,gamm,p,psi,beta,rg,rgol,rl,v,poc,scale;
   int type;

   printf("Please enter E(MeV):\n", &KE);
   scanf("%lf", &KE);
   printf("E=%lf\n",KE);

   printf("Please enter type:\n", &type);
   scanf("%d", &type);
   printf("type=%d\n",type);

   /* identify mass in kilogram */
   if(type==1) m = me; /* electron*/
   if(type==2) m = mp;
   else m = ma;
 /*  else m = 1.67262171e-27;*/ /* ion */

   E = KE + m;  /* MeV */
   gamm = E/m;  /* no unit */
   printf("gamm=%lf\n",gamm); 
   
   p = sqrt(E*E-m*m); /* MeV/c */
 /*  psi= p*(1.e6*1.60217653e-19/3.e8);*/ /* kg m/s */
   psi= p*(2.*q/300.);   /* alpha only */
   printf("p=%.10lf MeV/c = %le kg m/s\n",p,psi);

   beta = sqrt(1.-(1./(gamm*gamm))); /* v in unit of c */
   printf("beta=%.10lf\n",beta);
   rg = p/(2.*300.*B00*AU); /* AU */
   printf("rg=%.10lf\n",rg);
   rgol=rg/lperp;
   printf("rgol=%.10lf\n",rgol);      
   rl = (gamm*m*c*beta)/(q*B00*AU);
   printf("rl=%.10lf\n",rl);      

   v=beta*c;
   poc=(m*v)/(m*c);
   printf("poc=%.10lf\n",poc);
   
 /*  scale = (q*B0*lperp/gamm*psi);*/
   scale = (lperp*beta)/rg ; /* nondimension */
   printf("scale=%.15lf\n",scale);

}
