Platon Technologies
not logged in Login Registration
EnglishSlovak
open source software development celebrating 10 years of open source development! Saturday, October 17, 2026
[BACK] [TXT][DIR] Up to [Platon] / dead / diplomovka / xpdp / src

Diff for dead/diplomovka/xpdp/src/field.c between version 1.1 and 1.2

version 1.1, 2001/09/20 14:12:00 version 1.2, 2002/09/09 17:12:46
Line 32  void fields(void)
Line 32  void fields(void)
         // the numerical (implicit) susceptibility chi.  If          // the numerical (implicit) susceptibility chi.  If
         // explicit, sp_chi_scale=0 for all species => chi=0.          // explicit, sp_chi_scale=0 for all species => chi=0.
   
         for (j=0; j<=nc; j++) {          for (j = 0; j <= nc; j++) {
                 rho[j] = rhoback;                  rho[j] = rhoback;
                 chi[j] = 0.0;                  chi[j] = 0.0;
         }          }
   
         for (isp=0; isp<nsp; isp++) {          for (isp = 0; isp < nsp; isp++) {
                 for (j=0; j<=nc; j++) {                  for (j = 0; j <= nc; j++) {
                         chi[j] += sp_n[isp][j]*sp_chi_scale[isp];                          chi[j] += sp_n[isp][j] * sp_chi_scale[isp];
                         rho[j] += sp_n[isp][j]*q[isp]/dx;                          rho[j] += sp_n[isp][j] * q[isp] / dx;
                 }                  }
         } // }}}          } // }}}
         // recalculate the coeffs. for the tridiagnol field solve  {{{          // recalculate the coeffs. for the tridiagnol field solve  {{{
   
         for (j=1; j<nc; j++) {          for (j = 1; j < nc; j++) {
                 tri_a[j] = 1.0 + 0.5*(chi[j] + chi[j-1]);                  tri_a[j] = 1.0 + 0.5 * (chi[j] + chi[j - 1]);
                 tri_c[j] = 1.0 + 0.5*(chi[j] + chi[j+1]);                  tri_c[j] = 1.0 + 0.5 * (chi[j] + chi[j + 1]);
                 tri_b[j] = -(tri_a[j] + tri_c[j]);                  tri_b[j] = -(tri_a[j] + tri_c[j]);
         }          }
         if (bc_flag == BC1 || bc_flag == BC2) {          if (bc_flag == BC1 || bc_flag == BC2) {
Line 55  void fields(void)
Line 55  void fields(void)
                 tri_c[0] =  tri_a[1];                  tri_c[0] =  tri_a[1];
         }          }
         else if (bc_flag == BC4) {          else if (bc_flag == BC4) {
                 tri_b[0] = 2.25*extl/dt/dt + 1.5*extr/dt + 1/extc;                  tri_b[0] = 2.25 * extl / dt / dt + 1.5 * extr / dt + 1 / extc;
                 tri_b[0] = -(tri_a[1] + dx/tri_b[0]/area/epsilon);                  tri_b[0] = -(tri_a[1] + dx / tri_b[0] / area / epsilon);
                 tri_c[0] =  tri_a[1];                  tri_c[0] =  tri_a[1];
         } // }}}          } // }}}
   
         // solve a tridiagnal matrix to get Phi[i]  {{{          // solve a tridiagnal matrix to get Phi[i]  {{{
   
         jp = backj;          jp = backj;
         for (isp=0; isp<nsp; isp++) jp += jwall[isp];          for (isp = 0; isp < nsp; isp++)
                   jp += jwall[isp];
   
         if (nstrt) phi[0] = source();          if (nstrt)
                   phi[0] = source();
         bet = tri_b[nstrt];          bet = tri_b[nstrt];
         phi[nstrt] = (*bcptr)()/bet;          phi[nstrt] = (*bcptr)()/bet;
   
         for (j=nstrt+1; j<nc; j++) {          for (j = nstrt + 1; j < nc; j++) {
                 gam[j] = tri_c[j-1]/bet;                  gam[j] = tri_c[j-1]/bet;
                 bet = tri_b[j] - tri_a[j]*gam[j];                  bet = tri_b[j] - tri_a[j]*gam[j];
                 phi[j] = -(se*rho[j] + tri_a[j]*phi[j-1])/bet;                  phi[j] = -(se*rho[j] + tri_a[j]*phi[j-1])/bet;
         }          }
         phi[nc] = 0.;          phi[nc] = 0.;
         for (j=nc-2; j>= nstrt; j--) phi[j] -= gam[j+1]*phi[j+1];          for (j = nc - 2; j >= nstrt; j--)
                   phi[j] -= gam[j + 1] * phi[j + 1];
   
         // }}}          // }}}
   
         // calculate E(x) from phi(x)   {{{          // calculate E(x) from phi(x)   {{{
   
         s = 0.5/dx;          s = 0.5 / dx;
         for (j=1; j<nc; j++) e[j] = (phi[j-1] - phi[j+1])*s;          for (j = 1; j < nc; j++)
                   e[j] = (phi[j-1] - phi[j+1]) * s;
         e[0]  = (1.0+0.5*(chi[0]+chi[1]))*(phi[0] -phi[1])/dx  -rho[0]*0.5*dx/epsilon;  
         e[nc] = (1.0+0.5*(chi[nc-1]+chi[nc]))*(phi[nc-1] -phi[nc])/dx +rho[nc]*0.5*dx/epsilon;          e[0]  = (1.0 + 0.5 * (chi[0] + chi[1])) *
                   (phi[0] - phi[1]) / dx - rho[0] * 0.5 * dx / epsilon;
           e[nc] = (1.0 + 0.5 * (chi[nc-1] + chi[nc])) *
                   (phi[nc-1] - phi[nc]) / dx + rho[nc] * 0.5 * dx / epsilon;
         // }}}          // }}}
   
         // calculate the external circuit stuff          // calculate the external circuit stuff
Line 99  float source(void)
Line 105  float source(void)
 {  {
         if (dcramped) {          if (dcramped) {
                 if (t < risetime) {                  if (t < risetime) {
                         if (ramp) return (ramp*t);                          if (ramp)
                         else     return (.5*dcbias*(1- cos(t*w0)));                                  return (ramp*t);
                           else
                                   return (0.5 * dcbias *(1 - cos(t * w0)));
                 }                  }
                 else       return (dcbias);                  else
                           return (dcbias);
         }          }
         else return (dcbias + ramp*t + acbias*sin(t*w0 + theta0));          else
                   return (dcbias + ramp * t + acbias * sin(t * w0 + theta0));
 } // }}}  } // }}}
   
 // Fuction void field_init(void) {{{  // Fuction void field_init(void) {{{
Line 133  void field_init(void)
Line 143  void field_init(void)
                 circuitptr = circuit2;                  circuitptr = circuit2;
                 nstrt = 0;                  nstrt = 0;
         }          }
         else if (extl < 1e-30 && extr < 1e-30 && extc >= 1e5*epsilon*area/length) {          else if (extl < 1e-30 && extr < 1e-30 && extc >= 1e5 * epsilon * area / length) {
                 // When L=R=0.0 and C -> infinity, SHORT CIRCUIT                  // When L=R=0.0 and C -> infinity, SHORT CIRCUIT
                 bc_flag = BC3;                  bc_flag = BC3;
                 bcptr = bc3;                  bcptr = bc3;
Line 146  void field_init(void)
Line 156  void field_init(void)
                 circuitptr = circuit4;                  circuitptr = circuit4;
                 nstrt = 0;                  nstrt = 0;
         } // }}}          } // }}}
         se = dx*dx/epsilon;   // setting se = dx*dx/epsilon          se = dx * dx / epsilon;   // setting se = dx*dx/epsilon
 } // }}}  } // }}}
   
 // When C = 0: OPEN CIRCUIT  // When C = 0: OPEN CIRCUIT
Line 154  void field_init(void)
Line 164  void field_init(void)
 // Fuction float bc1(void) {{{  // Fuction float bc1(void) {{{
 float bc1(void)  float bc1(void)
 {  {
         return (-(sigma + dt*jp + 0.5*rho[0]*dx)*dx/epsilon);          return (-(sigma + dt * jp + 0.5 * rho[0] * dx) * dx / epsilon);
 } // }}}  } // }}}
 // Fuction void circuit1(void) {{{  // Fuction void circuit1(void) {{{
 void circuit1(void)  void circuit1(void)
Line 168  void circuit1(void)
Line 178  void circuit1(void)
 // Fuction float bc2(void) {{{  // Fuction float bc2(void) {{{
 float bc2(void)  float bc2(void)
 {  {
         return (-(sigma + dt*(source()/area + jp) + 0.5*rho[0]*dx)*dx/epsilon);          return (-(sigma + dt * (source() / area + jp) + 0.5 * rho[0] * dx) * dx / epsilon);
 } // }}}  } // }}}
 // Fuction void circuit2(void) {{{  // Fuction void circuit2(void) {{{
 void circuit2(void)  void circuit2(void)
 {  {
         exti    = source();          exti    = source();
         oldsigma= sigma;          oldsigma= sigma;
         sigma  += dt*(jp + exti/area);          sigma  += dt * (jp + exti / area);
 } // }}}  } // }}}
   
 // When R=L=0, C -> infinity: SHORT CIRCUIT  // When R=L=0, C -> infinity: SHORT CIRCUIT
Line 183  void circuit2(void)
Line 193  void circuit2(void)
 // Fuction float bc3(void) {{{  // Fuction float bc3(void) {{{
 float bc3(void)  float bc3(void)
 {  {
         return (-rho[1]*dx*dx/epsilon - tri_a[1]*source());          return (-rho[1] * dx * dx / epsilon - tri_a[1] * source());
 } // }}}  } // }}}
 // Fuction void circuit3() {{{  // Fuction void circuit3() {{{
 void circuit3()  void circuit3()
 {  {
         oldsigma = sigma;          oldsigma = sigma;
         sigma= epsilon*e[0];          sigma= epsilon*e[0];
         exti = area*((sigma - oldsigma)/dt- jp);          exti = area * ((sigma - oldsigma) / dt- jp);
 } // }}}  } // }}}
   
 // The General case  // The General case
Line 202  float bc4(void)
Line 212  float bc4(void)
         float k;          float k;
   
         if (!a0) {          if (!a0) {
                 a0 = 2.25*extl/dt/dt + 1.5*extr/dt + 1/extc;                  a0 =  2.25 * extl / dt / dt + 1.5 * extr / dt + 1 / extc;
                 a1 = -6*extl/dt/dt - 2*extr/dt;                  a1 = -6    * extl / dt / dt - 2   * extr / dt;
                 a2 = 5.5*extl/dt/dt + .5*extr/dt;                  a2 =  5.5  * extl / dt / dt + 0.5 * extr / dt;
                 a3 = -2*extl/dt/dt;                  a3 = -2    * extl / dt / dt;
                 a4 = .25*extl/dt/dt;                  a4 =  0.25 * extl / dt / dt;
         }          }
         k = (a1*extq + a2*extq_1 + a3*extq_2 + a4*extq_3)/a0;          k = (a1 * extq + a2 * extq_1 + a3 * extq_2 + a4 * extq_3) / a0;
         return (-(0.5*rho[0]*dx + sigma + dt*jp          return (-(0.5 * rho[0] * dx +
                                                                 + (source()/a0 - k - extq)/area)*dx/epsilon);                                  sigma + dt * jp +
                                   (source() / a0 - k - extq) / area) * dx / epsilon);
 } // }}}  } // }}}
 // Function void circuit4(void) {{{  // Function void circuit4(void) {{{
 void circuit4(void)  void circuit4(void)
Line 219  void circuit4(void)
Line 230  void circuit4(void)
         float k;          float k;
   
         if (!a0) {          if (!a0) {
                 a0 = 2.25*extl/dt/dt + 1.5*extr/dt + 1/extc;                  a0 =  2.25 * extl / dt / dt + 1.5 * extr / dt + 1 / extc;
                 a1 = -6*extl/dt/dt - 2*extr/dt;                  a1 = -6    * extl / dt / dt - 2   * extr / dt;
                 a2 = 5.5*extl/dt/dt + .5*extr/dt;                  a2 =  5.5  * extl / dt / dt + 0.5 * extr / dt;
                 a3 = -2*extl/dt/dt;                  a3 = -2    * extl / dt / dt;
                 a4 = .25*extl/dt/dt;                  a4 =  0.25 * extl / dt / dt;
         }          }
         k = (a1*extq + a2*extq_1 + a3*extq_2 + a4*extq_3)/a0;          k = (a1 * extq + a2 * extq_1 + a3 * extq_2 + a4 * extq_3) / a0;
         extq_3 = extq_2;          extq_3 = extq_2;
         extq_2 = extq_1;          extq_2 = extq_1;
         extq_1 = extq;          extq_1 = extq;
         extq = (source() - phi[0])/a0 - k;          extq   = (source() - phi[0]) / a0 - k;
   
         exti = (extq - extq_1)/dt;          exti = (extq - extq_1) / dt;
         oldsigma= sigma;          oldsigma = sigma;
         sigma  += (jp + exti/area)*dt;          sigma  += (jp + exti / area) * dt;
 } // }}}  } // }}}
   
 // Modeline for ViM {{{1  // Modeline for ViM {{{1
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
Platon Group <platon@platon.org> http://platon.org/
Copyright © 2002-2006 Platon Group
Site powered by Metafox CMS 		Go to Top