cktsens.c

linux平台下类似著名的电路板作图软件 Spice的源代码

				delta_iI[j] *= -1.0;
			}

#ifdef ASDEBUG
			DEBUG(2) {
				printf("Effect of negating matrix:\n");
				spPrint(delta_Y, 0, 1, 1);
				for (j = 0; j < size; j++)
					printf("%d: %g, %g\n", j,
						delta_I[j], delta_iI[j]);
			}
#endif

			/* XXX swap back to temp states ??   Naw ... */
			(void) sens_temp(sg, ckt);

#ifdef ASDEBUG
			DEBUG(1) {
				if (sens_getp(sg, ckt, &value)) {
					continue;
				}

				printf("New value in device: %g\n",
					value.rValue);
			}
#endif

			sens_load(sg, ckt, is_dc);

#ifdef ASDEBUG
			DEBUG(2) {
				printf("Effect of changing the parameter:\n");
				spPrint(delta_Y, 0, 1, 1);
				for (j = 0; j < size; j++)
					printf("%d: %g, %g\n", j,
						delta_I[j], delta_iI[j]);
			}
#endif
			/* Set the perturbed variable back to it's
			 * original value
			 */

			value.rValue = sg->value;
			sens_setp(sg, ckt, &value);
			(void) sens_temp(sg, ckt); /* XXX is this necessary? */

			/* Back to business . . . */

#ifdef ASDEBUG
			DEBUG(2)
				for (j = 0; j < size; j++)
					printf("          E [%d] = %20.15g\n",
						j, E[j]);
#endif

			/* delta_Y E */
			spMultiply(delta_Y, delta_I_delta_Y, E,
				delta_iI_delta_Y, iE);

#ifdef ASDEBUG
			DEBUG(2)
				for (j = 0; j < size; j++)
					printf("delta_Y * E [%d] = %20.15g\n",
						j, delta_I_delta_Y[j]);
#endif

			/* delta_I - delta_Y E */
			for (j = 0; j < size; j++) {
				delta_I[j] -= delta_I_delta_Y[j];
				delta_iI[j] -= delta_iI_delta_Y[j];
			}

#ifdef ASDEBUG
			DEBUG(2) {
				printf(">>> Y:\n");
				spPrint(Y, 0, 1, 1);
				for (j = 0; j < size; j++)
					printf("%d: %g, %g\n", j,
						delta_I[j], delta_iI[j]);
			}
#endif
			/* Solve; Y already factored */
			spSolve(Y, delta_I, delta_I, delta_iI, delta_iI);

#ifdef ASDEBUG
			DEBUG(2) {
				for (j = 1; j < size; j++) {

					if (sg->is_instparam)
						printf("%d/%s.%s = %g, %g\n",
							j,
							sg->instance->GENname,
						sg->ptable[sg->param].keyword,
						delta_I[j], delta_iI[j]);
					else
						printf("%d/%s:%s = %g, %g\n",
							j,
							sg->instance->GENname,
						sg->ptable[sg->param].keyword,
						delta_I[j], delta_iI[j]);

				}
			}
#endif

			/* delta_I is now equal to delta_E */

			if (is_dc) {
				if (sen_info->output_volt)
					output_values[n] = delta_I
						[sen_info->output_pos->number]
						- delta_I
						[sen_info->output_neg->number];
				else {
					output_values[n] = delta_I[branch_eq];
				}
				output_values[n] /= delta_var;
			} else {
				if (sen_info->output_volt) {
					output_cvalues[n].real = delta_I
						[sen_info->output_pos->number]
						- delta_I
						[sen_info->output_neg->number];
					output_cvalues[n].imag = delta_iI
						[sen_info->output_pos->number]
						- delta_iI
						[sen_info->output_neg->number];
				} else {
					output_cvalues[n].real =
						delta_I[branch_eq];
					output_cvalues[n].imag =
						delta_iI[branch_eq];
				}
				output_cvalues[n].real /= delta_var;
				output_cvalues[n].imag /= delta_var;
			}

			n += 1;

		}

		if (is_dc)
			nvalue.v.vec.rVec = output_values;
		else
			nvalue.v.vec.cVec = output_cvalues;

		value.rValue = freq;
		OUTpData(sen_data, &value, &nvalue);
		freq = inc_freq(freq, sen_info->step_type, step_size);

	}

	(*SPfrontEnd->OUTendPlot)((GENERIC *) sen_data);

	if (is_dc) {
		FREE(output_values);	/* XXX free various vectors */
	} else {
		FREE(output_cvalues);	/* XXX free various vectors */
	}

	spDestroy(delta_Y);
	FREE(delta_I);
	FREE(delta_iI);

	ckt->CKTrhs = E;
	ckt->CKTirhs = iE;
	ckt->CKTmatrix = Y;
	ckt->CKTbypass = bypass;

#ifdef notdef
	for (j = 0; j <= ckt->CKTmaxOrder + 1; j++) {
		if (ckt->CKTstates[j])
			FREE(ckt->CKTstates[j]);
		ckt->CKTstates[j] = save_states[j];
	}
#endif
#endif

	return OK;
}

double
inc_freq(freq, type, step_size)
	double	freq;
	int	type;
	double	step_size;
{
	if (type != LINEAR)
		freq *= step_size;
	else
		freq += step_size;

	return freq;
}

double
next_freq(type, freq, stepsize)
	int	type;
	double	freq, stepsize;
{
	double	s;

	switch (type) {
	case SENS_DC:
		s = 0;
		break;

	case SENS_LINEAR:
		s = freq + stepsize;
		break;

	case SENS_DECADE:
	case SENS_OCTAVE:
		s = freq * stepsize;
		break;
	}
	return s;
}

int
count_steps(type, low, high, steps, stepsize)
	int	type;
	double	low, high;
	int	steps;
	double	*stepsize;
{
	double	s;
	int	n;

	if (steps < 1)
		steps = 1;

	switch (type) {
	default:
	case SENS_DC:
		n = 0;
		s = 0;
		break;

	case SENS_LINEAR:
		n = steps;
		s = (high - low) / steps;
		break;

	case SENS_DECADE:
		if (low <= 0.0)
			low = 1e-3;
		if (high <= low)
			high = 10.0 * low;
		n = steps * log10(high/low) + 1.01;
		s = pow(10.0, 1.0 / steps);
		break;

	case SENS_OCTAVE:
		if (low <= 0.0)
			low = 1e-3;
		if (high <= low)
			high = 2.0 * low;
		n = steps * log(high/low) / M_LOG2E + 1.01;
		s = pow(2.0, 1.0 / steps);
		break;
	}

	if (n <= 0)
		n = 1;

	*stepsize = s;
	return n;
}

static int
sens_load(sg, ckt, is_dc)
	CKTcircuit	*ckt;
	sgen		*sg;
	int		is_dc;
{
	int	(*fn)( );

	error = 0;

	if (!is_dc)
		fn = DEVices[sg->dev]->DEVacLoad;
	else
		fn = DEVices[sg->dev]->DEVload;

	if (fn)
		error = (*fn)(sg->model, ckt);
	else
		return 1;

	return error;
}


static int
sens_temp(sg, ckt)
	CKTcircuit	*ckt;
	sgen		*sg;
{
	int	(*fn)( );

	error = 0;

	fn = DEVices[sg->dev]->DEVtemperature;

	if (fn)
		error = (*fn)(sg->model, ckt);
	else
		return 1;

	return error;
}

/* Get parameter value */
int
sens_getp(sg, ckt, val)
	CKTcircuit	*ckt;
	sgen		*sg;
	IFvalue		*val;
{
	int	(*fn)( );
	int	pid;

	error = 0;

	if (sg->is_instparam) {
		fn = DEVices[sg->dev]->DEVask;
		pid = DEVices[sg->dev]->DEVpublic.instanceParms[sg->param].id;
		if (fn)
			error = (*fn)(ckt, sg->instance, pid, val, NULL);
		else
			return 1;
	} else {
		fn = DEVices[sg->dev]->DEVmodAsk;
		pid = DEVices[sg->dev]->DEVpublic.modelParms[sg->param].id;
		if (fn)
			error = (*fn)(ckt, sg->model, pid, val, NULL);
		else
			return 1;
	}

	if (error) {
		if (sg->is_instparam)
			printf("GET ERROR: %s:%s:%s -> param %s (%d)\n",
				DEVices[sg->dev]->DEVpublic.name,
				sg->model->GENmodName,
				sg->instance->GENname,
				sg->ptable[sg->param].keyword, pid);
		else
			printf("GET ERROR: %s:%s:%s -> mparam %s (%d)\n",
				DEVices[sg->dev]->DEVpublic.name,
				sg->model->GENmodName,
				sg->instance->GENname,
				sg->ptable[sg->param].keyword, pid);
	}

	return error;
}

/* Get parameter value */
static int
sens_setp(sg, ckt, val)
	sgen		*sg;
	CKTcircuit	*ckt;
	IFvalue		*val;
{
	int	(*fn)( );
	int	pid;

	error = 0;

	if (sg->is_instparam) {
		fn = DEVices[sg->dev]->DEVparam;
		pid = DEVices[sg->dev]->DEVpublic.instanceParms[sg->param].id;
		if (fn)
			error = (*fn)(pid, val, sg->instance, NULL);
		else
			return 1;
	} else {
		fn = DEVices[sg->dev]->DEVmodParam;
		pid = DEVices[sg->dev]->DEVpublic.modelParms[sg->param].id;
		if (fn)
			error = (*fn)(pid, val, sg->model);
		else
			return 1;
	}

	if (error) {
		if (sg->is_instparam)
			printf("SET ERROR: %s:%s:%s -> param %s (%d)\n",
				DEVices[sg->dev]->DEVpublic.name,
				sg->model->GENmodName,
				sg->instance->GENname,
				sg->ptable[sg->param].keyword, pid);
		else
			printf("SET ERROR: %s:%s:%s -> mparam %s (%d)\n",
				DEVices[sg->dev]->DEVpublic.name,
				sg->model->GENmodName,
				sg->instance->GENname,
				sg->ptable[sg->param].keyword, pid);
	}

	return error;
}