nipzmuller.c

支持数字元件仿真的SPICE插件

/*
 * Copyright (c) 1985 Mani B. srivastava
 */

    /* NIpzMuller(ckt,ptrlistptr,type)
     *This is a routine implementing the Muller's method of 
     *finding roots of a complex polynomial
     *ptrlistptr is the pointer to a pointer of a list of roots
     *type is a number reflecting the type of analysis being done
     */

#include "prefix.h"
#include <stdio.h>
#include "PZdefs.h"
#include <math.h>
#include "util.h"
#include "complex.h"
#include "CKTdefs.h"
#include "DEVdefs.h"
#include "SPerror.h"
#include "suffix.h"

RCSID("NIpzMuller.c $Revision: 1.1 $ on $Date: 91/04/02 11:47:31 $")

int
NIpzMuller(ckt,ptrlistptr,type)
    register CKTcircuit *ckt;
    root **ptrlistptr;
    int type;
{
    int itr,error,NPASS;
    int power0,power1,power2,power3,power,scale;
    int oldscale,newscale;
    int scale1,scale2,scale0;
    int NIpzSolve();
    int first,num;
    double test;
    double temp1,temp2;
    double XMAX,PZTOL,PZABS,NMAX,NROOT,MATSIZE;
    double xmag,hmag,amax1,xi,xr;
    SPcomplex x0,x1,x2,x3,x;
    SPcomplex fx0,fx1,fx2,fx3,fx;
    SPcomplex h2,h1,h;
    SPcomplex lambda,lambda2;
    SPcomplex droot,denom1,denom;
    SPcomplex delta2;
    SPcomplex a,b,c;
    SPcomplex ctemp0,ctemp1,ctemp2,ctemp3,ctemp4,ctemp5;
    SPcomplex x2minusx,x1minusx,x0minusx;
    root *temproot;
    root *listptr;
    static char *itmsg = "Too many iterations without convergence";
    static char *magmsg = "Magnitude overflow in pole-zero analysis";
    /*
    printf("MULLER:entering \n");
    */
    listptr=(root *)(NULL);
    XMAX=1.0e20;
    PZTOL=1.0e-6;
    PZABS=1.0e-6;
    NPASS=1;
    NMAX=50;
    NROOT=0;
    MATSIZE=SMPmatSize(ckt->CKTmatrix);

    /*choose initial values x0 , x1 , x2 */
    lp100:
    first=0;
    scale=0;
    oldscale=0;
    x0.real= 1.763e-3;
    x0.imag=0.0;
    x1.real= -1.391e-5;
    x1.imag=0.0;
    x2.real= -1.271e0;
    x2.imag=0.0;
    /*calculate f(x0) , f(x1) , f(x2) 
     *note that f(s)=( det[Y(s)])/((s-r1)*...*(s-rk))
     *where r1,r2,r3...rk are the roots that have been found
     */
    error=NIpzSolve(ckt,&x2,listptr,&fx2,type,&power2);
    if (error) return(error);
    error=NIpzSolve(ckt,&x1,listptr,&fx1,type,&power1);
    if (error) return(error);
    error=NIpzSolve(ckt,&x0,listptr,&fx0,type,&power0);
    if (error) return(error);
    /*start search for a new root*/
    lp170:
    itr=3;
    lp175:
    /*first find a scaling factor for better numerical calculation*/
    if ((fabs(fx2.real)>0.0) || (fabs(fx2.imag)>0.0) ||
            (fabs(fx1.real)>0.0) || (fabs(fx2.imag)>0.0) ||
            (fabs(fx0.real)>0.0) || (fabs(fx2.imag)>0.0)) {
        newscale=0;
        num=0;
        if ((fabs(fx2.real)>0.0) || (fabs(fx2.imag)>0.0)) {
             /* printf("f(0)=(%5.6e,%5.6e) \n",fx2.real,fx2.imag); */
             scale2=1-power2;
             newscale=scale2;
             num++;
        }
        if ((fabs(fx1.real)>0.0) || (fabs(fx2.imag)>0.0)) {
             /* printf("f(0)=(%5.6e,%5.6e) \n",fx1.real,fx1.imag); */
             scale1=1-power1;
             num++;
             newscale=(newscale+scale1)/num;
        }
        if ((fabs(fx0.real)>0.0) || (fabs(fx2.imag)>0.0)) {
             /* printf("f(0)=(%5.6e,%5.6e) \n",fx0.real,fx0.imag); */
             scale0=1-power0;
             num++;
             newscale=((newscale*(num-1))+scale0)/num;
        }
     } else if (NPASS==1) {
        x3.real=1.0;
        x3.imag=0.0;
        AGAIN:
        error=NIpzSolve(ckt,&x3,listptr,&fx3,type,&power);
        if (error) return(error);
        /* printf("f(0)=(%5.6e,%5.6e) \n",fx3.real,fx3.imag); */
        if (fabs(fx3.real)>0.0) {
            newscale=1-power;
        } else {
            x3.real=x3.real+1.0;
            goto AGAIN;
        };
    };
    /* printf("newscale=%5.d \n",newscale); */
    scale=newscale-oldscale;
    oldscale=newscale;
    /* printf("scale=%5.d \n",scale); */
    if (first==0) {
        if ((fx2.real != 0.0) || (fx2.imag != 0.0)) {
            fx2.real=(pow(10.0,(double)(power2+scale))) * fx2.real;
            fx2.imag=(pow(10.0,(double)(power2+scale))) * fx2.imag;
        };
        if ((fx1.real != 0.0) || (fx1.imag != 0.0)) {
            fx1.real=(pow(10.0,(double)(power1+scale))) * fx1.real;
            fx1.imag=(pow(10.0,(double)(power1+scale))) * fx1.imag;
        };
        if ((fx0.real != 0.0) || (fx0.imag != 0.0)) {
            fx0.real=(pow(10.0,(double)(power0+scale))) * fx0.real;
            fx0.imag=(pow(10.0,(double)(power0+scale))) * fx0.imag;
        };
        first=1;
    } else {
        if ((fx2.real != 0.0) || (fx2.imag != 0.0)) {
            fx2.real=(pow(10.0,(double)(scale))) * fx2.real;
            fx2.imag=(pow(10.0,(double)(scale))) * fx2.imag;
        };
        if ((fx1.real != 0.0) || (fx1.imag != 0.0)) {
            fx1.real=(pow(10.0,(double)(scale))) * fx1.real;
            fx1.imag=(pow(10.0,(double)(scale))) * fx1.imag;
        };
        if ((fx0.real != 0.0) || (fx0.imag != 0.0)) {
            fx0.real=(pow(10.0,(double)(scale))) * fx0.real;
            fx0.imag=(pow(10.0,(double)(scale))) * fx0.imag;
        };
    };
    /* printf("x0=(%5.6e,%5.6e) \n",x0.real,x0.imag); */
    /* printf("fx0=(%5.6e,%5.6e) \n",fx0.real,fx0.imag); */
    /* printf("x1=(%5.6e,%5.6e) \n",x1.real,x1.imag); */
    /* printf("fx1=(%5.6e,%5.6e) \n",fx1.real,fx1.imag); */
    /* printf("x2=(%5.6e,%5.6e) \n",x2.real,x2.imag); */
    /* printf("fx2=(%5.6e,%5.6e) \n",fx2.real,fx2.imag); */

    /*now check whether f(x) is nearly a constant => all roots
     *have been found
     *in fact if this happens in any iteration , it means that
     *all roots have been found and we exit
     */
     /*calculate temp1=|f(x2)-f(x0)|*/
    temp1=fabs(fx2.real-fx0.real)+fabs(fx2.imag-fx0.imag);
    /* printf("temp1=%5.6e \n",temp1); */
    /*calculate temp2=|f(x2)-f(x1)|*/
    temp2=fabs(fx2.real-fx1.real)+fabs(fx2.imag-fx1.imag);
    /* printf("temp2=%5.6e\n",temp2); */
    /*check whether fx2,fx1,fx0 are nearly equal*/
    if ((temp1+temp2) > (PZTOL*fabs(fx2.real))) {
        goto lp200;
    };
    x3.real=100.0*x0.real;
    x3.imag=100.0*x0.imag;
    error=NIpzSolve(ckt,&x3,listptr,&fx3,type,&power3);
    if (error) return(error);
    if ((fx3.real != 0.0) || (fx3.imag != 0.0)) {
        fx3.real=(pow(10.0,(double)(power3+newscale))) * fx3.real;
        fx3.imag=(pow(10.0,(double)(power3+newscale))) * fx3.imag;
    };
    /*calculate temp1=|f(x2)-f(x3)|*/
    temp1=fabs(fx2.real-fx3.real)+fabs(fx2.imag-fx3.imag);
    /* printf("temp1a=%5.6e\n",temp1); */
    /* printf("temp2a=%5.6e\n",temp2); */
    if ((temp1+temp2) <= (PZTOL*fabs(fx2.real))) {
        goto lp1000;
    };
lp200:
    if (NROOT>MATSIZE) goto lp1000;
    /*calculate h2=x2-x1*/
    h2.real=x2.real - x1.real ;
    h2.imag=x2.imag - x1.imag ;
    /*calculate h1=x1-x0*/
    h1.real=x1.real - x0.real ;
    h1.imag=x1.imag - x0.imag ;
    /*calculate lambda2=h2 / h1*/
    lambda2.real = h2.real;
    lambda2.imag = h2.imag;
    DC_DIVEQ( &(lambda2.real) , &(lambda2.imag) , h1.real , h1.imag ) ;
    /*calculate delta2=1 + lambda2*/
    delta2.real=1.0+lambda2.real;
    delta2.imag=lambda2.imag;
    /*calculate ctemp0=lambda2*lambda2*/
    DC_MULT(lambda2.real , lambda2.imag , lambda2.real ,
        lambda2.imag , &(ctemp0.real) , &(ctemp0.imag));
    /*calculate ctemp1=fx0 * lambda2 * lambda2*/
    DC_MULT(fx0.real , fx0.imag , ctemp0.real , ctemp0.imag ,
        &(ctemp1.real) , &(ctemp1.imag));
    /*calculate ctemp2=fx1 * delta2*/
    DC_MULT(fx1.real , fx1.imag , delta2.real , delta2.imag ,
        &(ctemp2.real) , &(ctemp2.imag));
    /*calculate ctemp3=fx2 * lambda2*/
    DC_MULT(fx2.real , fx2.imag , lambda2.real , lambda2.imag ,
        &(ctemp3.real) , &(ctemp3.imag));
    /*calculate ctemp4=fx1 * delta2 * lambda2*/
    DC_MULT(ctemp2.real , ctemp2.imag , lambda2.real , lambda2.imag ,
        &(ctemp4.real) , &(ctemp4.imag));
    /*calculate a=fx0 * lambda2 ^ 2 
     *            - fx1 * lambda2 * delta2
     *            + fx2 * lambda2
     */
    a.real=ctemp1.real+ctemp3.real-ctemp4.real;
    a.imag=ctemp1.imag+ctemp3.imag-ctemp4.imag;
    /*calculate c=fx2 * delta2*/
    c.real=fx2.real+ctemp3.real;
    c.imag=fx2.imag+ctemp3.imag;
    /*calculate b=fx0 * lambda2 ^ 2 
     *            - fx1 * delta2 ^ 2