radtr.c

一维等离子体静电粒子模拟程序

float lorentz(float x)
{
  return(1/(1+x*x));
}
float doppler(float x)
{
  return(exp(-x*x));
}
float voigt(float x)		/* NOT PREPARED YET */
{
  return(1/(1+x*x));
}


    /********** The integral omega for Diffusion Coeff. ***********/
float omega(float x)
{
  int   i;
  float tstar[27] = {0.250, 0.267, 0.286, 0.308, 0.333, 0.364, 0.400, 0.444,
                     0.500, 0.571, 0.667, 0.800, 1.000, 1.330, 2.000, 2.500,
                     2.860, 3.330, 4.000, 5.000, 6.670, 10.00, 12.50, 16.67,
                     25.00, 50.00, 100.0},
        omega[27] = {1.8714, 1.8344, 1.7962, 1.7560, 1.7142, 1.6702, 1.6238,
                     1.5748, 1.5228, 1.4676, 1.4090, 1.3466, 1.2802, 1.2096,
                     1.1348, 1.1042, 1.0888, 1.0738, 1.0592, 1.0450, 1.0318,
                     1.0198, 1.0154, 1.0112, 1.0072, 1.0034, 1.000};
  float temp1, temp2=0;

  if      (x<tstar[0])  temp2=omega[0];
  else if (x>tstar[26]) temp2=omega[26];
  else
    for (i=1;i<27;i++)
      if (x<tstar[i]) {
        temp1=(x-tstar[i-1])/(tstar[i]-tstar[i-1]);
        temp2=omega[i]*temp1+omega[i-1]*(1-temp1);
      }
  return(temp2);
}

/*******************************************************************/
void calculate_Akm()
{
  register int i, j, k;
  float x, k_shape, ri, rju, rjl, temp;

  for (i=0;i<7;i++) {
    I1_l[i] /=pow(3.75,2*i);
    I0_l[i]/=pow(3.75,2*i);
  }
  for (i=0;i<9;i++) {
    I1_u[i] *=pow(3.75,i);
    I0_u[i]*=pow(3.75,i);
  }
  for (i=0;i<7;i++) {
    I0_l[i] /=pow(3.75,2*i);
    I1_l[i]/=pow(3.75,2*i);
  }
  for (i=0;i<7;i++) {
    K0_l[i] /=pow(2.0,2*i);
    K1_l[i]/=pow(2.0,2*i);
  }
  for (i=0;i<7;i++) {
    K0_u[i] *=pow(2.0,i);
    K1_u[i]*=pow(2.0,i);
  }
  /*------------------ frequency integration --------------------*/
  for (i=0;i<nc_RT;i++) {
    for (j=0;j<nc_RT;j++) {
      A[i][j]=0;
      for (k=0;k<=(xbin-1)/2;k++) {
        x=x_array[k];   
        k_shape=k0*k_function(x)*Dr;
        ri=(i+0.5)*k_shape;
        rju=(j+1.0)*k_shape;
        rjl=j*k_shape;
        
        if (i>j) {
          if (rjl<3.75) temp=(phi_integral1(i,j+1,k_shape)-
                              phi_integral1(i,j,k_shape))/M_PI;
          else          temp=bK0I1(ri,rju) - bK0I1(ri,rjl);
        }
        else if (i<j) {
          if (ri <3.75) temp=(phi_integral2(i,j+1,k_shape)-
                              phi_integral2(i,j,k_shape))/M_PI;
          else          temp=aK1I0(rjl,ri) - aK1I0(rju,ri);
        }
        else {
          if (rjl<3.75) temp=(phi_integral2(i,j+1,k_shape)-
                              phi_integral1(i,j,k_shape))/M_PI;
          else          temp= bK0I1(ri,ri) - bK0I1(ri,rjl) 
                            + aK1I0(ri,ri) - aK1I0(rju,ri);
        }
        if (k==0 || k==(xbin-1)/2) temp*=0.5;
        A[i][j]+=temp*k_function(x);
      }                      
      A[i][j]*=2*Cx*dfreqx;
/*      if (fabs(A[i][j])<1e-10) A[i][j]=0;*/
    }
  } 

}

/*******************************************************************/
float bK0I1(float a, float b)
{
  float result, z, h1, temp, temp2;
  register int i, j;

  z=a-b;
  for (i=0;i<15;i++) K0uI1u[i]=0;
  temp=1;
  for (i=0;i<7;i++) {
    for (j=0;j<9;j++) K0uI1u[i+j]+=K0_u[i]*I1_u[j]*temp;
    temp*=b/a;                                          /* temp=pow(b/a,i) */
  }
  
  {
    /*------- if alpha=beta -------*/
    if (fabs(z)<1e-15) {			
      h1=0.0;
      temp2=b;
      for (i=0;i<15;i++) {
        h1+=K0uI1u[i]/(i+1.0)/temp2;
        temp2*=b;
      }
      result=h1;
    }
    /*------- if alpha>beta -------*/
    else {				
      h1=0.0;
      temp=exp(-z)-z*ei1m(z);			/* when i=0, temp=Ei2(z) */
      temp2=b;
      for (i=0;i<15;i++) {
        h1+=K0uI1u[i]*temp/temp2;
        temp=(exp(-z)-z*temp)/(i+2.0);	/* Ei_j+1(z)=1/j*(exp(-z)-z*Ei_j(z)) */
        temp2*=b;
      }
      result=h1/sqrt(a/b);
    }
  }  

  return (result*b);
}

/*******************************************************************/
float aK1I0(float a, float b)
{
  float result, z, h1, temp, temp2;
  register int i,j;

  z=a-b;
  for (i=0;i<15;i++) K1uI0u[i]=0;
  temp=1;
  for (i=0;i<7;i++) {
    for (j=0;j<9;j++) K1uI0u[i+j]+=K1_u[i]*I0_u[j]*temp;
    temp*=b/a;                                          /* temp=pow(b/a,i) */
  }
  {
    /*------- if alpha=beta -------*/
    if (fabs(z)<1e-15) {     
      h1=0.0;
      temp2=b;
      for (i=0;i<15;i++) {
        h1+=K1uI0u[i]/(i+1.0)/temp2;
        temp2*=b;
      }
      result=h1;
    }
    /*------- if alpha>beta -------*/
    else {
      h1=0.0;
      temp=exp(-z)-z*ei1m(z);                   /* when i=0, temp=Ei2(z) */
      temp2=b;
      for (i=0;i<15;i++) {
        h1+=K1uI0u[i]*temp/temp2;
        temp=(exp(-z)-z*temp)/(i+2.0);  /* Ei_j+1(z)=1/j*(exp(-z)-z*Ei_j(z)) */
        temp2*=b;
      }
      result=h1/sqrt(a/b);
    }
  }
  return (result*a);
}

/*------------------- With the method I, for k>=m -------------------*/
float phi_integral1(int k, int m, float k0kx)    
{
  register int i;
  float theta_max, dtheta, theta, ha, hb, result=0;
  int   Nmax=256;
 
  if (m==0) return(0);

  theta_max=asin(m/r_grid2[k]);
  Nmax*=theta_max/M_PI;
  if (Nmax<64) Nmax=64;
  dtheta=theta_max/(float)Nmax;
  for (i=1;i<Nmax;i++) {
     theta=dtheta*i;
     ha=r_grid2[k]*cos(theta)+sqrt(m*m-pow(r_grid2[k]*sin(theta),2));
     hb=r_grid2[k]*cos(theta)-sqrt(m*m-pow(r_grid2[k]*sin(theta),2));
     result+= Ki2(k0kx*hb)-Ki2(k0kx*ha);
  }  
  ha=r_grid2[k]+m;
  hb=r_grid2[k]-m;
  result += 0.5*(Ki2(k0kx*hb)-Ki2(k0kx*ha));
  result *= dtheta;
/*  printf("C %e %f %f %e %e\n",k0kx,ha,hb,dtheta,result); */
  return(result);
}

/*------------------- With the method I, for k<m --------------------*/
float phi_integral2(int k, int m, float k0kx)    
{
  register int i;
  float theta_max, dtheta, theta, ha, result=0;
  int   Nmax=256;
  
  theta_max=M_PI;
  dtheta=theta_max/(float)Nmax;
  for (i=1;i<Nmax;i++) {
     theta=dtheta*i;
     ha=r_grid2[k]*cos(theta)+sqrt(m*m-pow(r_grid2[k]*sin(theta),2));
     result+= 1-Ki2(k0kx*ha);
  }
  ha=r_grid2[k]+m;
  result+= 0.5*(1-Ki2(k0kx*ha));
  ha=-r_grid2[k]+m;
  result+= 0.5*(1-Ki2(k0kx*ha));
  result*= dtheta;
/*  printf("%e %e\n",k0kx,result); */
  return(result);
}


/*******************************************************************/

float ei1m(float z)
{
  float a0,a1,a2,a3,a4,a5, b1,b2,b3,b4, ei1;

/*  if (z<=0) {
    puts("Ei1(z) function cannot be evaluated for z<=0");
    exit(-1);
  }*/

  if (z<1) {
    a0=-0.57721566; a1=0.99999193;
    a2=-0.24991055; a3=0.05519968;
    a4=-0.00976004; a5=0.00107857;
    ei1=a0+z*(a1+z*(a2+z*(a3+z*(a4+z*a5)))) - log(z);
  }
  else {
    a1=8.5733287401; a2=18.0590169730;
    a3=8.6347608925; a4=0.2677737343;
    b1=9.5733223454; b2=25.6329561486;
    b3=21.0996530827; b4=3.9584969228;
    ei1=exp(-z)/z*(a4+z*(a3+z*(a2+z*(a1+z))))/(b4+z*(b3+z*(b2+z*(b1+z))));
  }
  return (ei1);
}

float GE0(float z) 			/* GE0(z) = sqrt(PI)*erfc(sqrt(z)) */
{
  float t, result;
/*---------------------- epsilon < 1.5e-7 ----------------------*
  t=1.0/(1+0.3275911*sqrt(z));
  result = t*(0.254829592 +t*(-0.284496736 + t*(1.421413741 + t*(-1.453152027 
         + t*1.061405429))))*exp(-z);
*/
/*---------------------- epsilon < 2.5e-5 ----------------------*/
  t=1.0/(1+0.47047*sqrt(z));
  result = t*(0.3480242 + t*(-0.0958798 + t*0.7478556))*exp(-z);

  return (sqrt(M_PI)*result); 
}
      
/* for GNU Plot
Error Function for x>0, 
plot 1.0/(1+0.3275911*x)*(0.254829592 +1.0/(1+0.3275911*x)*(-0.284496736 + 1.0/(1+0.3275911*x)*(1.421413741 + 1.0/(1+0.3275911*x)*(-1.453152027 + 1.0/(1+0.3275911*x)*1.061405429))))*exp(-x*x), 1.0/(1+0.47047*x)*(0.3480242 + 1.0/(1+0.47047*x)*(-0.0958798 + 1.0/(1+0.47047*x)*0.7478556))*exp(-x*x)
*/

float J0(float x)
{
  float y, theta0, result;
  if (x<3) {
    y=x/3;
    y*=y;
    result= 1.0 + y*(-2.2499997 + y*(1.2656208 + y*(-0.3163866 + y*(0.0444479 +
                  y*(-0.0039444 + y*0.0002100)))));
  }
  else {
    y=3./x;
    theta0= x - 0.78539816 + y*(-0.04166397 + y*(-0.00003954 + y*(0.00262573 +
            y*(-0.00054125 + y*(-0.00029333 + y*0.00013558)))));
    result= cos(theta0)*(0.79788456 + y*(-0.00000077 + y*(-0.00552740 +
            y*(-0.00009512 +y*(0.00137237 + y*(-0.00072805 + y*0.00014476))))))
            /sqrt(x);
  }
  return (result);
}

/***************************************************************/
float Ki2(float z)
{
  float result;

  if (z<1e-9) result = 1.0;
  else if (z>30) result=0;
  else         result = zK1(z)-z*exp(-z)*Ki1(z);
  return (result);
}

float Ki1(float z)				/* exp(z) Ki1(z) */
/* This function computes the second repeated integral of
   the modified Bessel function second kind, zero order
   by means of a recursion formula (see Abramowitz)      */
{
  float result, t, u, h1, tsq;

  if (z<1e-12) result=M_PI/2;
  else if (z<1.0) {
    t=z/2;
    tsq=t*t;
    u=log(t)+0.5772156649;
    h1=(1-u) + tsq*( (4.0/9.0-u/3.0) + tsq*( (0.085-0.05*u) +
       tsq*( (0.007842025-u/252.0) + tsq*(4.23311036e-4 - u/5184.0) )));
    result=(-2*t*h1+M_PI/2)*exp(z);
  }
  else if (z<3.8) {
    t=1/z;
    result=1.0/sqrt(z)*(1.23711091 +t*(-0.62015283 + t*(0.44927294 +
           t*(-0.23121965 + t*(0.11001445 + t*(-0.09732933 + t*0.04438948))))));
  }
  else {
    t=7.0/z;
    result=1.0/sqrt(z)*(1.25331414 + t*(-0.11190289 + t*(0.02576646 +
           t*(-0.00933994 + t*(0.00417454 + t*(-0.00163271 + t*0.00033934))))));
  }
  return (result);
}

float zK1(float z)                              /* exp(z) Ki-1(z) */
{
  float t, result;

  if (z<2) {
    t=z/2.0;
    t*=t;
    result= z*log(z/2.0)*I1(z) + 1 + t*(0.15443144 + t*(-0.67278579 +
        t*(-0.18156897 + t*(-0.01919402 + t*(-0.00110404 - t*0.00004686)))));
  }
  else {
    t=2.0/z;
    result=exp(-z)*sqrt(z)*( 1.25331414 + t*(0.23498619 + t*(-0.03655620 +
           t*(0.01504268 + t*(-0.00780353 + t*(0.00325614 - t*0.00068245) )))));
  }
  return (result);
}

float I1(float z)
{
  float t, result;
  if (z<=3.75) {
    t=z/3.75;
    t*=t;
    result= z*(0.5 + t*(0.87890594 + t*(0.51498869 + t*(0.15084934 +
               t*(0.02658733 + t*(0.00301532 + t*0.00032411))))));
  }
  else {
    t=3.75/z;
    result= exp(z)/sqrt(z)*(0.39894228 + t*(-0.03988024 + t*(-0.00362018 +
           t*(0.00163801 +t*(-0.01031555 + t*(0.02282967 + t*(-0.02895312 +
           t*(0.01787654 - t*0.00420059))))))));
  }
  return (result);
}