cplload.c

ngspice又一个电子CAD仿真软件代码.功能更全

      for (k = 0; k < noL; k++)       /**  node voltage  **/
     for (l = 0; l < noL; l++)    /**  different mode  **/
        if (cp->h3t[j][k][l]) {
           if (cp->h3t[j][k][l]->ifImg) {
          double er, ei, a, b, a1, b1, a2, b2;
          TMS *tms;
          tms = cp->h3t[j][k][l];
          expC(tms->tm[1].x, tms->tm[2].x, h, &er, &ei);
          a2 = h1 * tms->tm[1].c;
          b2 = h1 * tms->tm[2].c;
          a = tms->tm[1].cnv_i;
          b = tms->tm[2].cnv_i;
           multC(a, b, er, ei, &a, &b);
          multC(a2, b2, (double) v1_i[l][k] * er + v2_i[l][k], (double) v1_i[l][k] * ei, &a1, &b1);
          tms->tm[1].cnv_i = a + a1;
          tms->tm[2].cnv_i = b + b1;
          a = tms->tm[1].cnv_o;
          b = tms->tm[2].cnv_o;
          multC(a, b, er, ei, &a, &b);
          multC(a2, b2, (double) v1_o[l][k] * er + v2_o[l][k], (double) v1_o[l][k] * ei, &a1, &b1);
          tms->tm[1].cnv_o = a + a1;
          tms->tm[2].cnv_o = b + b1;
          tm = &(tms->tm[0]);
          e = exp((double) tm->x * h);
          tm->cnv_i = tm->cnv_i * e + h1 * tm->c * 
                  (v1_i[l][k] * e + v2_i[l][k]);
          tm->cnv_o = tm->cnv_o * e + h1 * tm->c * 
                  (v1_o[l][k] * e + v2_o[l][k]);
          ff[j] += tms->aten * v2_o[l][k] + tm->cnv_o +
             2.0 * tms->tm[1].cnv_o;
             gg[j] += tms->aten * v2_i[l][k] + tm->cnv_i +
            2.0 * tms->tm[1].cnv_i;
           } else {
          for (i = 0; i < 3; i++) {   /**     3 poles     **/
             tm = &(cp->h3t[j][k][l]->tm[i]);
             e =  exp((double) tm->x * h);
             tm->cnv_i = tm->cnv_i * e + h1 * tm->c * (v1_i[l][k] * e + v2_i[l][k]);
             tm->cnv_o = tm->cnv_o * e + h1 * tm->c * (v1_o[l][k] * e + v2_o[l][k]);
             ff[j] += tm->cnv_o;
             gg[j] += tm->cnv_i;
          }

          ff[j] += cp->h3t[j][k][l]->aten * v2_o[l][k];
          gg[j] += cp->h3t[j][k][l]->aten * v2_i[l][k];
           }
        }
      for (k = 0; k < noL; k++)       /**  line current  **/
     for (l = 0; l < noL; l++)    /**  different mode  **/
        if (cp->h2t[j][k][l]) {
           if (cp->h2t[j][k][l]->ifImg) {
          double er, ei, a, b, a1, b1, a2, b2;
          TMS *tms;
          tms = cp->h2t[j][k][l];
          expC(tms->tm[1].x, tms->tm[2].x, h, &er, &ei);
          a2 = h1 * tms->tm[1].c;
          b2 = h1 * tms->tm[2].c;
          a = tms->tm[1].cnv_i;
          b = tms->tm[2].cnv_i;
          multC(a, b, er, ei, &a, &b);
          multC(a2, b2, (double) i1_i[l][k] * er + i2_i[l][k], (double) i1_i[l][k] * ei, &a1, &b1);
          tms->tm[1].cnv_i = a + a1;
          tms->tm[2].cnv_i = b + b1;
          a = tms->tm[1].cnv_o;
          b = tms->tm[2].cnv_o;
          multC(a, b, er, ei, &a, &b);
          multC(a2, b2, (double) i1_o[l][k] * er + i2_o[l][k], (double) i1_o[l][k] * ei, &a1, &b1);
          tms->tm[1].cnv_o = a + a1;
          tms->tm[2].cnv_o = b + b1;
          tm = &(tms->tm[0]);
          e = exp((double) tm->x * h);
          tm->cnv_i = tm->cnv_i * e + h1 * tm->c * 
                  (i1_i[l][k] * e + i2_i[l][k]);
          tm->cnv_o = tm->cnv_o * e + h1 * tm->c * 
                  (i1_o[l][k] * e + i2_o[l][k]);
          ff[j] += tms->aten * i2_o[l][k] + tm->cnv_o +
             2.0 * tms->tm[1].cnv_o;
             gg[j] += tms->aten * i2_i[l][k] + tm->cnv_i +
            2.0 * tms->tm[1].cnv_i;
           } else {
          for (i = 0; i < 3; i++) {   /**     3 poles     **/
             tm = &(cp->h2t[j][k][l]->tm[i]);
             e =  exp((double) tm->x * h);
             tm->cnv_i = tm->cnv_i * e + h1 * tm->c * (i1_i[l][k] * e + i2_i[l][k]);
             tm->cnv_o = tm->cnv_o * e + h1 * tm->c * (i1_o[l][k] * e + i2_o[l][k]);
             ff[j] += tm->cnv_o;
             gg[j] += tm->cnv_i;
          }

          ff[j] += cp->h2t[j][k][l]->aten * i2_o[l][k];
          gg[j] += cp->h2t[j][k][l]->aten * i2_i[l][k];
           }
        }
   } 

   for (i = 0; i < noL; i++) {
	  *(ckt->CKTrhs + l1[i]) = ff[i];
	  *(ckt->CKTrhs + l2[i]) = gg[i];
   }

   return(ext);
}

static int 
update_cnv(CPLine *cp, float h)
{
   int i, j, k;
   int noL;
   double ai, bi, ao, bo;
   double e, t;
   TMS *tms;
   TERM *tm;

   noL = cp->noL;
   for (j = 0; j < noL; j++)
      for (k = 0; k < noL; k++) {
     ai = cp->in_node[k]->V;
     ao = cp->out_node[k]->V;
     bi = cp->in_node[k]->dv;
     bo = cp->out_node[k]->dv;

     if (cp->h1t[j][k]) {
        if (cp->h1t[j][k]->ifImg) {
           tms = cp->h1t[j][k];
           if (tms == NULL)
          continue;
           tm = &(tms->tm[0]);
           e = cp->h1e[j][k][0];
           t = tm->c / tm->x;
           update_cnv_a(tms, h, ai, ao, ai - bi * h, ao - bo * h,
              cp->h1e[j][k][1], cp->h1e[j][k][2]);
           bi *= t;
           bo *= t;
           tm->cnv_i = (tm->cnv_i - bi*h) * e + (e - 1.0)*(ai*t +
              1.0e+12*bi/tm->x);
           tm->cnv_o = (tm->cnv_o - bo*h) * e + (e - 1.0)*(ao*t + 
              1.0e+12*bo/tm->x);
        } else 
           for (i = 0; i < 3; i++) {
          tm = &(cp->h1t[j][k]->tm[i]);

          e = cp->h1e[j][k][i];
    
          t = tm->c / tm->x;
          bi *= t;
          bo *= t;

          tm->cnv_i = (tm->cnv_i - bi*h) * e + (e - 1.0)*(ai*t + 1.0e+12*bi/tm->x);
          tm->cnv_o = (tm->cnv_o - bo*h) * e + (e - 1.0)*(ao*t + 1.0e+12*bo/tm->x);
           }  
     }
      }
   return 0;
}

static VI_list
*new_vi(void)
{
   VI_list *q;

	if (pool_vi) {
		q = pool_vi;
  	    pool_vi = pool_vi->pool;
		return(q);
    } else return((VI_list *) malloc (sizeof (VI_list)));
}

static void
free_vi(VI_list *q)
{
  	q->pool = pool_vi;
	pool_vi = q;
}


static int 
add_new_vi(CPLinstance *here, CKTcircuit *ckt, int time)
{
   VI_list *vi;
   int i, noL;
   CPLine *cp, *cp2;

   cp = here->cplines;
   cp2 = here->cplines2;

   vi = new_vi();
   vi->time = time;
   noL = cp->noL;
   for (i = 0; i < noL; i++) {
	  /* 
      vi->v_i[i] = cp->in_node[i]->V;
      vi->v_o[i] = cp->out_node[i]->V;
	  */
      vi->v_i[i] = *(ckt->CKTrhsOld + here->CPLposNodes[i]);
      vi->v_o[i] = *(ckt->CKTrhsOld + here->CPLnegNodes[i]);
      vi->i_i[i] = *(ckt->CKTrhsOld + here->CPLibr1[i]);
      vi->i_o[i] = *(ckt->CKTrhsOld + here->CPLibr2[i]);
   }
   cp->vi_tail->next = vi;
   cp2->vi_tail->next = vi;
   vi->next = NULL;
   cp->vi_tail = vi;
   cp2->vi_tail = vi;

   return(1);
}


static int
get_pvs_vi(t1, t2, cp, v1_i, v2_i, i1_i, i2_i, v1_o, v2_o, i1_o, i2_o)
   CPLine *cp;
   int t1, t2;
   double v1_i[MAX_CP_TX_LINES][MAX_CP_TX_LINES];
   double v2_i[MAX_CP_TX_LINES][MAX_CP_TX_LINES];
   double i1_i[MAX_CP_TX_LINES][MAX_CP_TX_LINES];
   double i2_i[MAX_CP_TX_LINES][MAX_CP_TX_LINES]; 
   double v1_o[MAX_CP_TX_LINES][MAX_CP_TX_LINES]; 
   double v2_o[MAX_CP_TX_LINES][MAX_CP_TX_LINES];
   double i1_o[MAX_CP_TX_LINES][MAX_CP_TX_LINES]; 
   double i2_o[MAX_CP_TX_LINES][MAX_CP_TX_LINES];
{
   double ta[MAX_CP_TX_LINES], tb[MAX_CP_TX_LINES]; 
   VI_list *vi, *vi1;
   double f;
   int i, j;
   int  mini = -1;
   double minta = 123456789.0;
   int ext = 0;
   int noL;

   noL = cp->noL;

   for (i = 0; i < noL; i++) {
      ta[i] = t1 - cp->taul[i];
      tb[i] = t2 - cp->taul[i];
      if (ta[i] < minta) {
		 minta = ta[i];
		 mini = i;
      }
   }

   for (i = 0; i < noL; i++) {

      ratio[i] = 0.0;

      if (tb[i] <= 0) {
		 for (j = 0; j < noL; j++) {
			i1_i[i][j] = i2_i[i][j] = i1_o[i][j] = i2_o[i][j] = 0.0;
			v1_i[i][j] = v2_i[i][j] = cp->dc1[j];
			v1_o[i][j] = v2_o[i][j] = cp->dc2[j];
		 }
      } else {
		 if (ta[i] <= 0) {
			for (j = 0; j < noL; j++) {
			   i1_i[i][j] = i1_o[i][j] = 0.0; 
			   v1_i[i][j] = cp->dc1[j];
			   v1_o[i][j] = cp->dc2[j]; 
			}
			vi1 = cp->vi_head;
			vi = vi1->next;
		 } else {
			vi1 = cp->vi_head;
			for (vi = vi1->next; vi->time < ta[i]; ) {
			   /* if (i == mini)
			  		 free_vi(vi1); */
			   vi1 = vi;

			   /* new */ 
			   vi = vi->next;
			   if (vi == NULL) goto errordetect;
			}
			f = (ta[i] - vi1->time) / (vi->time - vi1->time);
			for (j = 0; j < noL; j++) {
			   v1_i[i][j] = vi1->v_i[j] + f * (vi->v_i[j] - vi1->v_i[j]);
			   v1_o[i][j] = vi1->v_o[j] + f * (vi->v_o[j] - vi1->v_o[j]);
			   i1_i[i][j] = vi1->i_i[j] + f * (vi->i_i[j] - vi1->i_i[j]);
			   i1_o[i][j] = vi1->i_o[j] + f * (vi->i_o[j] - vi1->i_o[j]);
			}
			if (i == mini)
			   cp->vi_head = vi1;
		 }

		 if (tb[i] > t1) {
			/*
			fprintf(stderr, "pvs: time = %d\n", t2);
			 */
			ext = 1;

			ratio[i] = f = (tb[i] - t1) / (t2 - t1);

			if (vi)
			   for (; vi->next; vi = vi->next);
			else
			   vi = vi1;
			f = 1 - f;
			for (j = 0; j < noL; j++) {
			   v2_i[i][j] = vi->v_i[j] * f;
			   v2_o[i][j] = vi->v_o[j] * f;
			   i2_i[i][j] = vi->i_i[j] * f;
           	   i2_o[i][j] = vi->i_o[j] * f;
            }
     	  } else {
			for (; vi->time < tb[i];) {
			   vi1 = vi;

			   /* new */
			   vi = vi->next;
			   if (vi == NULL) goto errordetect;
			} 
		  
			f = (tb[i] - vi1->time) / (vi->time - vi1->time);
			for (j = 0; j < noL; j++) {
			   v2_i[i][j] = vi1->v_i[j] + f * (vi->v_i[j] - vi1->v_i[j]);
			   v2_o[i][j] = vi1->v_o[j] + f * (vi->v_o[j] - vi1->v_o[j]);
			   i2_i[i][j] = vi1->i_i[j] + f * (vi->i_i[j] - vi1->i_i[j]);
			   i2_o[i][j] = vi1->i_o[j] + f * (vi->i_o[j] - vi1->i_o[j]);
			}
     	}
      }
   }

   return(ext);

errordetect:
	fprintf(stderr, "your maximum time step is too large for tau.\n");
	fprintf(stderr, "decrease max time step in .tran card and try again\n");
	exit(0);
}


static int 
update_delayed_cnv(CPLine *cp, float h)
{
   int i, j, k;
   float *ratio;
   double f;
   VI_list *vi;
   TMS *tms;
   int noL;

   h *= 0.5e-12;
   ratio = cp->ratio;
   vi = cp->vi_tail;
   noL = cp->noL;

   for (k = 0; k < noL; k++)     /*  mode  */
      if (ratio[k] > 0.0) 
     for (i = 0; i < noL; i++)  /*  current eqn  */
        for (j = 0; j < noL; j++) {
           tms = cp->h3t[i][j][k];
           if (tms == NULL)
          continue;
           f = h * ratio[k] * vi->v_i[j];
           tms->tm[0].cnv_i += f *  tms->tm[0].c;
           tms->tm[1].cnv_i += f *  tms->tm[1].c;
           tms->tm[2].cnv_i += f *  tms->tm[2].c;

           f = h * ratio[k] * vi->v_o[j];
           tms->tm[0].cnv_o += f *  tms->tm[0].c;
           tms->tm[1].cnv_o += f *  tms->tm[1].c;
           tms->tm[2].cnv_o += f *  tms->tm[2].c;

           tms = cp->h2t[i][j][k];
           f = h * ratio[k] * vi->i_i[j];
           tms->tm[0].cnv_i += f *  tms->tm[0].c;
           tms->tm[1].cnv_i += f *  tms->tm[1].c;
           tms->tm[2].cnv_i += f *  tms->tm[2].c;

           f = h * ratio[k] * vi->i_o[j];
           tms->tm[0].cnv_o += f *  tms->tm[0].c;
           tms->tm[1].cnv_o += f *  tms->tm[1].c;
           tms->tm[2].cnv_o += f *  tms->tm[2].c;
        }
   return(1);
}


static int 
expC(double ar, double ai, double h, double *cr, double *ci)
{
	double e, cs, si;

	e = exp((double) ar * h);
	cs = cos((double) ai * h);
	si = sin((double) ai * h);
	*cr = e * cs;
	*ci = e * si;

	return(1);
}

static int 
multC(double ar, double ai, double br, double bi, 
      double *cr, double *ci)
{
	double tp;

	tp = ar*br - ai*bi;
	*ci = ar*bi+ai*br;
	*cr = tp;

	return (1);

}

static void
update_cnv_a(TMS *tms, float h, double ai, double ao, double bi, double bo, 
             double er, double ei)
{
   double a, b, a1, b1;

   h *= 0.5e-12;
   multC(tms->tm[1].c, tms->tm[2].c, er, ei, &a1, &b1);
   multC(tms->tm[1].cnv_i, tms->tm[2].cnv_i, er, ei, &a, &b);
   tms->tm[1].cnv_i = a + h * (a1 * bi + ai * tms->tm[1].c);
   tms->tm[2].cnv_i = b + h * (b1 * bi + ai * tms->tm[2].c);

   multC(tms->tm[1].cnv_o, tms->tm[2].cnv_o, er, ei, &a, &b);
   tms->tm[1].cnv_o = a + h * (a1 * bo + ao * tms->tm[1].c);
   tms->tm[2].cnv_o = b + h * (b1 * bo + ao * tms->tm[2].c);
}

static int 
divC(double ar, double ai, double br, double bi, double *cr, double *ci)
{
	double t;

	t = br*br + bi*bi;
	*cr = (ar*br + ai*bi) / t;
	*ci = (ai*br - ar*bi) / t;

	return(1);
}