JSOC_Documentation: proj/globalhs/libs/pkbgn/fixrts.c Source File

00001 // $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/pkbgn/fixrts.c,v 1.1 2013/04/28 07:50:36 tplarson Exp $
00002 
00003 #include <math.h>
00004 #include "complex.h"
00005 
00006 #define NPMAX 100
00007 #define ZERO Complex(0.0,0.0)
00008 #define ONE Complex(1.0,0.0)
00009 
00010 void fixrts(d,npoles)
00011 float d[];
00012 int npoles;
00013 {
00014     int i,j,polish;
00015     fcomplex a[NPMAX],roots[NPMAX];
00016     void zroots();
00017 
00018     a[npoles]=ONE;
00019     for (j=npoles-1;j>=0;j--)
00020         a[j]=Complex(-d[npoles-j],0.0);
00021     polish=1;
00022     zroots(a,npoles,roots,polish);
00023     for (j=1;j<=npoles;j++)
00024         if (Cabs(roots[j]) > 1.0)
00025             roots[j]=Cdiv(ONE,Conjg(roots[j]));
00026     a[0]=Csub(ZERO,roots[1]);
00027     a[1]=ONE;
00028     for (j=2;j<=npoles;j++) {
00029         a[j]=ONE;
00030         for (i=j;i>=2;i--)
00031             a[i-1]=Csub(a[i-2],Cmul(roots[j],a[i-1]));
00032         a[0]=Csub(ZERO,Cmul(roots[j],a[0]));
00033     }
00034     for (j=0;j<=npoles-1;j++)
00035         d[npoles-j] = -a[j].r;
00036 }