nipzmuller.c

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

     *            + fx2 * lambda2
     *            + fx2 * delta2
     */
    b.real=ctemp1.real-ctemp2.real+ctemp3.real-ctemp4.real+c.real;
    b.imag=ctemp1.imag-ctemp2.imag+ctemp3.imag-ctemp4.imag+c.imag;
    /*calculate ctemp5=b^2 - 4 * a * c */
    ctemp5.real=b.real * b.real - b.imag * b.imag -
        4.0 * a.real * c.real + 4.0 * a.imag * c.imag;
    ctemp5.imag=2.0 * b.real * b.imag -
        4.0 * a.real * c.imag - 4.0 * a.imag * c.real;
/*now find the complex square root of ctemp5*/
    if (ctemp5.real == 0.0) {
        if (ctemp5.imag == 0.0) {
            droot.real = 0.0;
            droot.imag = 0.0;
        } else if (ctemp5.imag > 0.0) {
            droot.real = sqrt (0.5 * ctemp5.imag);
            droot.imag = droot.real;
        } else { 
            droot.imag = sqrt( -0.5 * ctemp5.imag);
            droot.real= - droot.imag;
        }
    } else if (ctemp5.real > 0.0) {
        if (ctemp5.imag == 0.0) {
            droot.real = sqrt(ctemp5.real);
            droot.imag = 0.0;
        } else if (ctemp5.imag < 0.0) {
            droot.real = -sqrt(0.5 * (sqrt(ctemp5.real * ctemp5.real
                    + ctemp5.imag * ctemp5.imag) + ctemp5.real));
        } else {
            droot.real =  sqrt(0.5 * (sqrt(ctemp5.real * ctemp5.real
                + ctemp5.imag * ctemp5.imag) + ctemp5.real));
        }
        droot.imag = ctemp5.imag / (2.0 * droot.real);
    } else { /* ctemp5.real < 0.0) */
        if (ctemp5.imag == 0.0) {
            droot.real = 0.0;
            droot.imag = sqrt(- ctemp5.real);
        } else {
            if (ctemp5.imag < 0.0) {
                droot.imag = -sqrt(0.5 * (sqrt(ctemp5.real * ctemp5.real
                        + ctemp5.imag * ctemp5.imag) - ctemp5.real));
            } else {
                droot.imag =  sqrt(0.5 * (sqrt(ctemp5.real * ctemp5.real
                        + ctemp5.imag * ctemp5.imag) - ctemp5.real));
                droot.real = ctemp5.imag / (2.0 * droot.imag);
            }
        }
    };
    /*calculate denom=b+sqrt(b^2-4ac) */
    denom.real=b.real + droot.real;
    denom.imag=b.imag + droot.imag;
    /*calculate denom1=b-sqrt(b^2-4ac) */
    denom1.real=b.real - droot.real;
    denom1.imag=b.imag - droot.imag;
    if ((fabs(denom1.real) + fabs(denom1.imag)) > 
            (fabs(denom.real) + fabs(denom.imag))) {
        denom.real = denom1.real;
        denom.imag = denom1.imag;
    };
    /*calculate smallest zero lambda = -2c/denom */
    lambda.real= -2.0 * c.real;
    lambda.imag= -2.0 * c.imag;
    DC_DIVEQ(&(lambda.real) , &(lambda.imag) , denom.real , denom.imag);
    /*calculate h=lambda * h2 */
    DC_MULT(lambda.real , lambda.imag , h2.real , h2.imag ,&(h.real),&(h.imag));
    /*calculate new x=x2 + h */
    x.real=x2.real + h.real;
    x.imag=x2.imag + h.imag;

    /*do convergence test*/
    xmag=fabs(x.real) + fabs(x.imag);
    if (xmag > XMAX) {
        if (NPASS >= 2) {
            errMsg = MALLOC(strlen(magmsg)+1);
            strcpy(errMsg,magmsg);
            return(E_MAGEXCEEDED);
        } else {
            NPASS++;
            goto lp100;
        };
    };
    hmag=fabs(h.real) + fabs(h.imag);
    amax1=(xmag>PZABS) ? xmag : PZABS ;
    if (hmag <= PZTOL * amax1) goto lp300 ;
    error=NIpzSolve(ckt,&x,listptr,&fx,type,&power);
    if (error) return(error);
    if ((fx.real != 0.0) || (fx.imag != 0.0)) {
        fx.real=(pow(10.0,(double)(power+newscale))) * fx.real;
        fx.imag=(pow(10.0,(double)(power+newscale))) * fx.imag;
    } else {
        /*fx=(0,0)*/
        power=0;
    };
    /* printf("x=(%5.6e,%5.6e) \n",x.real,x.imag); */
    /* printf("fx=(%5.6e,%5.6e) \n",fx.real,fx.imag); */
    if ((fabs(fx.real) + fabs(fx.imag)) < 1.0e-25) {
        goto lp300 ; /*converged as |f(x2)|<1.0e-30 so that x2
               is a root*/
    };

    /*no convergence*/
    if (itr >= NMAX) {
        errMsg = MALLOC(strlen(itmsg)+1);
        strcpy(errMsg,itmsg);
        return(E_ITERLIM) ;
    }
    itr +=1 ;
    fx0.real = fx1.real ;
    fx0.imag = fx1.imag ;
    x0.real  = x1.real  ;
    x0.imag  = x1.imag  ;
    power0=power1;
    fx1.real = fx2.real ;
    fx1.imag = fx2.imag ;
    x1.real  = x2.real  ;
    x1.imag  = x2.imag  ;
    power1=power2;
    fx2.real = fx.real ;
    fx2.imag = fx.imag ;
    x2.real  = x.real  ;
    x2.imag  = x.imag  ;
    power2=power;
    goto lp175 ;

    /*converged*/
lp300:
    xr=x.real ;
    xi=x.imag ;
    xi=(fabs(xi) < (PZTOL * fabs(xr))) ? 0.0 : xi ;
    xi=(fabs(xi) < 1.0e-15) ? 0.0 : xi ;
    xr=(fabs(xr) < (PZTOL * fabs(xi))) ? 0.0 : xr ;
    xr=(fabs(xr) < 1.0e-15) ? 0.0 : xr ;
    /* printf("root=(%5.6e,%5.6e) \n",xr,xi); */
    NROOT++ ;
    /*store the root*/
    temproot=(root *) malloc((unsigned) sizeof(root));
    if (temproot==(root *)(NULL)) return(E_NOMEM);
    temproot->real=xr;
    temproot->imag=xi;
    temproot->next=listptr;
    listptr=temproot;
    x2minusx.real=x2.real-x.real;
    x2minusx.imag=x2.imag-x.imag;
    x1minusx.real=x1.real-x.real;
    x1minusx.imag=x1.imag-x.imag;
    x0minusx.real=x0.real-x.real;
    x0minusx.imag=x0.imag-x.imag;
    DC_DIVEQ(&(fx2.real),&(fx2.imag),x2minusx.real,x2minusx.imag);
    DC_DIVEQ(&(fx1.real),&(fx1.imag),x1minusx.real,x1minusx.imag);
    DC_DIVEQ(&(fx0.real),&(fx0.imag),x0minusx.real,x0minusx.imag);
    if (xi==0.0) goto lp400;
    /*store the conjugate root*/
    NROOT++ ;
    temproot=(root *) malloc((unsigned) sizeof(root));
    if (temproot==(root *)(NULL)) return(E_NOMEM);
    temproot->real=xr;
    temproot->imag= -xi;
    /* printf("root=(%5.6e,%5.6e) \n",xr,-xi); */
    temproot->next=listptr;
    listptr=temproot;
    x2minusx.real=x2.real-x.real;
    x2minusx.imag=x2.imag+x.imag;
    x1minusx.real=x1.real-x.real;
    x1minusx.imag=x1.imag+x.imag;
    x0minusx.real=x0.real-x.real;
    x0minusx.imag=x0.imag+x.imag;
    DC_DIVEQ(&(fx2.real),&(fx2.imag),x2minusx.real,x2minusx.imag);
    DC_DIVEQ(&(fx1.real),&(fx1.imag),x1minusx.real,x1minusx.imag);
    DC_DIVEQ(&(fx0.real),&(fx0.imag),x0minusx.real,x0minusx.imag);
    lp400:
    /*normalise fx2*/
    test=(fabs(fx2.real)>fabs(fx2.imag)) ?
    fabs(fx2.real) : fabs(fx2.imag) ;
    num=0;
    if (test !=0.0) {
        while (test>1.0) {
            test=test/10.0;
            fx2.real=fx2.real/10.0;
            fx2.imag=fx2.imag/10.0;
            num=num + 1;
        };
        while (test<=0.1) {
            test=test*10.0;
            fx2.real=fx2.real*10.0;
            fx2.imag=fx2.imag*10.0;
            num=num - 1;
        };
    };
    power2=num-oldscale;
    /*normalise fx1*/
    test=(fabs(fx1.real)>fabs(fx1.imag)) ?
    fabs(fx1.real) : fabs(fx1.imag) ;
    num=0;
    if (test !=0.0) {
        while (test>1.0) {
            test=test/10.0;
            fx1.real=fx1.real/10.0;
            fx1.imag=fx1.imag/10.0;
            num=num + 1;
        };
        while (test<=0.1) {
            test=test*10.0;
            fx1.real=fx1.real*10.0;
            fx1.imag=fx1.imag*10.0;
            num=num - 1;
        };
    };
    power1=num-oldscale;
    /*normalise fx0*/
    test=(fabs(fx0.real)>fabs(fx0.imag)) ?
    fabs(fx0.real) : fabs(fx0.imag) ;
    num=0;
    if (test !=0.0) {
        while (test>1.0) {
            test=test/10.0;
            fx0.real=fx0.real/10.0;
            fx0.imag=fx0.imag/10.0;
            num=num + 1;
        };
        while (test<=0.1) {
            test=test*10.0;
            fx0.real=fx0.real*10.0;
            fx0.imag=fx0.imag*10.0;
            num=num - 1;
        };
    };
    power0=num-oldscale;
    first=0;
    oldscale=0;
    goto lp170;
    
lp1000:
    *ptrlistptr=listptr;
    return(OK);
}