一个简单的CAD仿真程序——《电子电路的计算机辅助分析与设计方法(第2版)》附录源代码(略有修改)
注意，这个代码只支持简单的直流分析。
测试代码如下：
V1 1 0 5 R1 4 0 10 R2 3 0 100 R3 4 3 2 D1 1 2 0.01 1e-14 0.6 R4 3 2 10 z
头文件如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define R_Leq 10e-1 /*resistor value of inductor*/
#define R_Ceq 1.0e+7 /*resistor value of capacitor*/
void read_net_list();
void parser_net_list();
void LU_decomposition();
void output_net_list();
void output_operating_points();
void calculate_resistors(int);
void calculate_voltage_sources(int);
void calculate_current_sources(int);
void calculate_vccs(int);
void calculate_diode_operating_point(int);
int converge_check();
void calcluate_one_diode(int);
int get_max_mode(void);
static int
NUM_NODE,
NUM_DIM,
NUM_UNIT,
NUM_LINEAR,
NUM_CS,
NUM_NEW,
NUM_DIODE;
struct VOLTAGE_NODE
{
char name[5];
int node_i;
int node_j;
double value;
} NODES[10];
struct CS_NODE
{
char name[5];
int node_i;
int node_j;
int node_k;
int node_l;
double value;
}CS[10];
struct DIODE_NODE
{
char name[5];
int node_i;
int node_j;
double rs;
double is;
double vd;
double long id;
double long gd;
int line_num;
}DIODES[10];
struct CURRENT_NODE
{
int line_num;
int node_i;
int node_j;
}I_NODE[10];
源代码如下：
#include "CAD.h"
double A[30][30],B[30],LB[30];
int global_iter = 0, gloable_itermax = 150;
double VT = 0.026, ERR = 0.0001;
FILE *fp_input;
void initial();
void calculate_diodes();
/************************************************
主函数，处理文件输入或者屏幕输入，
以z或者Z作为文件结尾。
元件符号名不能超过4个字符，
目前只支持R、C、L、D、I、V和G这个几种简单的符号，
且只支持计算直流静态工作点，
即电感近似短路，电容近似开路。
作为一个简单的示例，最多支持30个节点（包括地）。
代码改编自《电子电路的计算机辅助分析与设计方法(第2版)》，
杨华中，罗嵘，汪慧编著；北京：清华大学出版社，2008.2
************************************************/
int main(int argc, char* argv[])
{
int return_value = 0;
if(argc == 2)
{
fp_input = fopen(argv[1], "r");
if(fp_input == NULL)
{
printf("Cannot open file: %s\n", argv[1]);
exit(0);
}
}
else
{
printf("Please input data by screen\n");
fp_input = stdin;
}
read_net_list(); // 读入网表
do
{
initial();// 初始化参数
parser_net_list(); // 解析网表
calculate_diodes(); // 计算二极管的静态工作点，获取电压和电流参数
LU_decomposition(); // LU分解计算AX=B的矩阵方程
return_value = converge_check(); // 检查迭代是否收敛
if(NUM_DIODE == 0)
return_value = 0;
}
while(return_value);
output_net_list(); // 输出一次读入的网表
output_operating_points(); // 输出求解的静态工作点
system("pause");
return 0;
}
/***************************************************************
导纳矩阵是A
输出电流矩阵是B
节点电压是X
矩阵方程是A*X=B
采用LU分解，即A=LU，来计算X，也就是节点电压法计算。
***************************************************************/
/************************************************
初始化A和B数组
************************************************/
void initial()
{
int i, j;
for(i = 0; i < 30; i++)
{
for (j = 0; j < 30; j++)
{
A[i][j] =0.0;
B[i] = 0.0;
}
}
NUM_NEW = 0;
}
/************************************************
循环处理所有的二极管参数
************************************************/
void calculate_diodes()
{
int m = 0;
while(m < NUM_DIODE)
{
calcluate_one_diode(m);
m++;
}
}
/************************************************
从文件或者屏幕读入电路的网表文件
************************************************/
void read_net_list()
{
int i=0, j=0, l=0, flag=0;
char symbol_name[5];
double var;
while(flag == 0)
{
fscanf(fp_input,"%s", symbol_name);
NUM_UNIT++;
switch(symbol_name[0])
{
default:
{
printf("\n There is no such component\n");
break;
}
case 'Z':
case 'z':
{
flag = 1;
NUM_UNIT--;
break;
}
case 'L':case 'l':
case 'C':case 'c':
{
fscanf(fp_input, "%d%d",&(NODES[i].node_i), &(NODES[i].node_j));
strcpy(NODES[i].name, symbol_name);
i++;
break;
}
case 'R':case 'r':
case 'I':case 'i':
case 'V':case 'v':
{
fscanf(fp_input, "%d%d%lf", &(NODES[i].node_i), &(NODES[i].node_j), &var);
strcpy(NODES[i].name, symbol_name);
NODES[i].value = var;
i++;
break;
}
case 'G':
case 'g':
{
fscanf(fp_input,"%d%d%d%d%lf",&(CS[j].node_k), &(CS[j].node_l), &(CS[j].node_i), &(CS[j].node_j), &var);
strcpy(CS[j].name, symbol_name);
CS[j].value = var;
j++;
break;
}
case 'D':
case 'd':
{
DIODES[l].is=1e-14;
DIODES[l].rs=1.0;
DIODES[l].vd=0.0;
fscanf(fp_input, "%d%d%lf%lf%lf", &(DIODES[l].node_i), &(DIODES[l].node_j), &(DIODES[l].rs),
&(DIODES[l].is), &(DIODES[l].vd));
strcpy(DIODES[l].name, symbol_name);
if(DIODES[l].is > 1e9) DIODES[l].is = 1e-14;
l++;
break;
}
}
}
NUM_LINEAR =i;
NUM_CS = j;
NUM_DIODE = l;
for(i =0 ; i < NUM_DIODE; i++)
{
DIODES[i].node_i--;
DIODES[i].node_j--;
}
for(i = 0; i < NUM_LINEAR; i++)
{
NODES[i].node_i--;
NODES[i].node_j--;
}
for(i = 0; i < NUM_CS; i++)
{
CS[i].node_i--;
CS[i].node_j--;
CS[i].node_k--;
CS[i].node_l--;
}
}
/************************************************
按照符号的首字母解析网表文件
************************************************/
void parser_net_list()
{
NUM_NODE = get_max_mode();
NUM_DIM = NUM_NODE;
for(int m = 0; m < NUM_LINEAR; m++)
{
switch(NODES[m].name[0])
{
case 'R':case 'r':
case 'C':case 'c':
{
calculate_resistors(m);
break;
}
case 'I':case 'i':
{
calculate_current_sources(m);
break;
}
case 'V':case 'v':
{
calculate_voltage_sources(m);
break;
}
}
}
for(int m = 0; m < NUM_CS; m++)
calculate_vccs(m);
}
/************************************************
获取电路的最大节点数目
************************************************/
int get_max_mode()
{
int n =0;
for(int m=0; m < NUM_LINEAR; m++)
{
if(n < NODES[m].node_i) n = NODES[m].node_i;
if(n < NODES[m].node_j) n = NODES[m].node_j;
}
for(int m = 0; m < NUM_CS; m++)
{
if( n < CS[m].node_i) n = CS[m].node_i;
if( n < CS[m].node_j) n = CS[m].node_j;
if( n < CS[m].node_k) n = CS[m].node_k;
if( n < CS[m].node_l) n = CS[m].node_l;
}
for(int m=0; m < NUM_DIODE; m++)
{
if(n < DIODES[m].node_i) n = DIODES[m].node_i;
if(n < DIODES[m].node_j) n = DIODES[m].node_j;
}
return n;
}
/************************************************
计算电阻的参数，将电阻变为电导，
注意本节点计算时，电导为正；
相邻节点计算时，电导为负。
列出电导的矩阵。
************************************************/
void calculate_resistors(int m)
{
int x, y;
double var;
if(NODES[m].name[0] == 'L'||NODES[m].name[0] == 'l') NODES[m].value = R_Leq;
if(NODES[m].name[0] == 'C'||NODES[m].name[0] == 'C') NODES[m].value = R_Ceq;
x = NODES[m].node_i;
y = NODES[m].node_j;
if(NODES[m].value >= 1.0e-9) var = 1.0/NODES[m].value;
else
{
printf("Component %s is short-cut\n", NODES[m].name);
exit(0);
}
if(x != -1)
{
A[x][x] = A[x][x] + var;
if( y != -1)
{
A[y][x] = A[y][x] - var;
A[x][y] = A[x][y] - var;
A[y][y] = A[y][y] + var;
}
}
else
{
if(y != -1)
A[y][y] = A[y][y] + var;
}
}
/************************************************
计算电流源的情况
************************************************/
void calculate_current_sources(int m)
{
if(NODES[m].node_i != -1)
B[NODES[m].node_i] = -NODES[m].value;
if(NODES[m].node_j != -1)
B[NODES[m].node_j] = +NODES[m].value;
}
/************************************************
计算电压源的情况
************************************************/
void calculate_voltage_sources(int m)
{
int x = NODES[m].node_i;
int y = NODES[m].node_j;
int flag = 0;
/* if at the same two nodes there is already a voltage source and the new added voltage source's value is not the same as the old one,
then the circuit data has error. Because the same tow nodes cannot has different voltage*/
for(int z = 0; z <= NUM_NEW-1; z++)
{
if(I_NODE[z].node_i == x && I_NODE[z].node_j == y
|| I_NODE[z].node_i == y && I_NODE[z].node_j == x)
{
if(fabs(B[I_NODE[z].line_num]-NODES[m].value)>1.0e-3)
{
printf("Error in circuit design on V\n");
exit(0);
}
else
{
flag = 1;
break;
}
}
}
if(flag == 0)
{
NUM_DIM++;
I_NODE[NUM_NEW].line_num = NUM_DIM;
I_NODE[NUM_NEW].node_i = x;
I_NODE[NUM_NEW].node_j = y;
NUM_NEW++;
B[NUM_DIM] = B[NUM_DIM] + NODES[m].value;
if(x != -1)
{
A[x][NUM_DIM] = A[x][NUM_DIM] + 1;
A[NUM_DIM][x] = A[NUM_DIM][x] + 1;
}
if(y != -1)
{
A[y][NUM_DIM] = A[y][NUM_DIM] - 1;
A[NUM_DIM][y] = A[NUM_DIM][y] - 1;
}
}
}
/************************************************
计算电压控制电压源的情况
************************************************/
void calculate_vccs(int m)
{
int i, j, k ,l;
double var;
i = CS[m].node_i;
j = CS[m].node_j;
k = CS[m].node_k;
l = CS[m].node_l;
var = CS[m].value;
if(i != -1 && k != -1)
A[k][i] = A[k][i] + var;
if(i != -1 && l != -1)
A[l][i] = A[l][i] - var;
if(j != -1 && k != -1)
A[k][j] = A[k][j] - var;
if(j != -1 && l != -1)
A[l][j] = A[l][j] + var;
}
/************************************************
计算二极管的情况
************************************************/
void calcluate_one_diode(int l)
{
calculate_diode_operating_point(l);
NUM_DIM++;
DIODES[l].line_num = NUM_DIM;
int y = DIODES[l].line_num;
double gd = DIODES[l].gd;
double id =DIODES[l].id;
double rs = DIODES[l].rs;
int il = DIODES[l].node_i;
int jl = DIODES[l].node_j;
if(il != -1)
{
A[il][il] = A[il][il] + 1/rs;
A[il][y] -= 1/rs;
A[y][il] -= 1/rs;
}
A[y][y] += 1/rs + gd;
if(jl != -1)
{
A[y][jl] -= gd;
A[jl][y] -= gd;
A[jl][jl] += gd;
B[jl] += id;
}
B[y] -= id;
}
/************************************************
计算二极管的电流和电压
************************************************/
void calculate_diode_operating_point(int l)
{
double diq = 1e-5;
double dis = DIODES[l].is;
double vq = log(1+diq/dis)*VT;
int i = DIODES[l].line_num;
double vd = LB[i];
int j = DIODES[l].node_j;
if(j != -1) vd -= LB[j];
if(global_iter == 0) vd = DIODES[l].vd;
if(global_iter == 0 && vd == 0.0) vd = vq;
if(global_iter != 0)
{
double id = DIODES[l].gd * vd + DIODES[l].id;
if(vd > vq && id >= diq)
vd = VT * log(id/dis+1);
else
if(vd > vq) vd = vq;
}
double gd = 0.0;
if(vd / VT > -100) gd = dis/VT * exp(vd/VT);
if(vd/VT > -100 && gd < 1e-14) gd = 1e-14;
double dj = gd*(VT-vd)-dis;
DIODES[l].gd = gd;
DIODES[l].id = dj;
DIODES[l].vd = vd;
}
/************************************************
LU方法分解矩阵求解
************************************************/
void LU_decomposition()
{
int Line_num[10], Col_num[10];
int MAX_line, MAX_col, NR;
int i, j, k, x;
double var, sum, MAX_a;
NR = NUM_DIM;
for( i = 0; i <= NR; i++)
{
Line_num[i] = i;
Col_num[i] = i;
}
for(k = 0; k <= NR-1; k++)
{
MAX_a = A[k][k];
MAX_line = k;
MAX_col = k;
for(i = k; i <= NR; i++)
{
for(j = k; j <= NR; j++)
{
var = fabs(A[i][i]);
if(var > MAX_a)
{
MAX_a = var;
MAX_line = i;
MAX_col = j;
}
}
}
for(j = 0; j <= NR; j++)
{
var = A[k][j];
A[k][j]= A[MAX_line][j];
A[MAX_line][j] = var;
}
x = Line_num[k];
Line_num[k] = Line_num[MAX_line];
Line_num[MAX_line] = x;
for( i = 0; i <= NR; i++)
{
var = A[i][k];
A[i][k] = A[i][MAX_col];
A[i][MAX_col] =var;
}
x = Col_num[k];
Col_num[k] = Col_num[MAX_col];
Col_num[MAX_col] = x;
for(i = k +1; i <= NR; i++)
A[i][k] = A[i][k]/A[k][k];
for(i = k+1; i <= NR; i++)
{
for(j =k+1; j <=NR; j++)
A[i][j] = A[i][j] - A[i][k]*A[k][j];
}
}
for(i = 0; i <= NR-1; i++)
{
x=Line_num[i];
var = B[i];
B[i] = B[x];
B[x] = var;
for(j = i+1; j <= NR; j++)
if(Line_num[j] == i) Line_num[j] =x;
}
for(i = 1; i<= NR; i++)
{
sum = 0;
for(j = 0; j <= i-1; j++)
sum = sum + A[i][j]*B[j];
B[i] = B[i] -sum;
}
if(A[NR][NR] < 1.e-6 && A[NR][NR] > -1.0e-6)
{
printf("Error in result: no exact answer\n");
exit(0);
}
B[NR]=B[NR]/A[NR][NR];
for(i = 1; i <= NR; i++)
{
sum = 0;
for(j = NR-i+1; j <= NR; j++)
sum = sum + A[NR-i][j]*B[j];
if(A[NR-i][NR-i] < 1e-6 && A[NR-i][NR-i] > -1e-3)
{
printf("Error in result: no exact answer\n");
exit(0);
}
B[NR-i] = (B[NR-i]-sum)/A[NR-i][NR-i];
}
for(i =0; i <= NR; i++)
Line_num[i] = Col_num[i];
for(i = 0; i <= NR; i++)
{
for(j = i; j <= NR; j++)
{
if(Line_num[j] == i)
{
var = B[i];
B[i] = B[j];
B[j] = var;
Line_num[j] = Line_num[i];
}
}
}
}
/************************************************
检查计算结果的收敛性
************************************************/
int converge_check()
{
int i;
double vk;
double ebsr = 0.0;
for( i = 0; i < 30; i++)
{
vk = fabs(LB[i] - B[i]);
if(ebsr < vk) ebsr = vk;
LB[i] = B[i];
}
global_iter++;
if(ebsr <= ERR)
{
printf("result: NR total iter=\t %d\n", global_iter);
return 0;
}
if(ebsr > ERR)
{
if(global_iter > gloable_itermax)
{
printf("Iter>%d, not converge\n", gloable_itermax);
return 0;
}
}
return 1;
}
/************************************************
输出网表，确认读入数据是否正确
************************************************/
void output_net_list()
{
printf("The total number of component is: %d\n", NUM_UNIT);
printf("the number of line2 is :%d\n", NUM_LINEAR);
for(int i = 0; i < NUM_LINEAR;i++)
printf("%s\t%d\t%d\t%f\n",NODES[i].name,NODES[i].node_i+1, NODES[i].node_j+1, NODES[i].value);
printf("the number of cs is: %d\n", NUM_CS);
for(int i = 0; i < NUM_CS; i++)
printf("%s\t%d\t%d\t%d\t%d\t%f\n",CS[i].name, CS[i].node_i+1, CS[i].node_j+1,
CS[i].node_k+1, CS[i].node_l+1, CS[i].value);
for(int i = 0; i < NUM_DIODE; i++)
printf("%s\t%d\t%d\t%f\t%f\n", DIODES[i].name, DIODES[i].node_i+1, DIODES[i].node_j+1, DIODES[i].rs, DIODES[i].is);
}
/************************************************
输出节点电压和节点之间的电流
************************************************/
void output_operating_points()
{
for(int i = 0; i <= NUM_NODE; i++)
printf("V%d = %g\n", i+1, B[i]);
for(int i = 0; i < NUM_NEW; i++)
printf("I(%d,%d)=%g\n",I_NODE[i].node_i+1, I_NODE[i].node_j+1, B[NUM_NODE+i+1]);
}
代码测试运行结果：
使用LTspice的仿真结果：
