break;
}
case3:
case4:
{
n2=2*m;
n1=n2+1;
work=malloc(n1*n1*sizeof(double));
w2=tan(pi*fhn);
w=w2-w1;
w0=w1*w2;
if(band==4)
{
for(i=0;i<=m/2;i++)
{
tmp=d[i];
d[i]=d[m-i];
d[m-i]=tmp;
tmp=c[i];
c[i]=c[m-i];
c[m-i]=tmp;
}
}
for(i=0;i<=n2;i++)
{
work[0*n1+i]=0.0;
work[1*n1+i]=0.0;
}
for(i=0;i<=m;i++)
{
tmpd=d[i]*pow(w,(m-i));
tmpc=c[i]*pow(w,(m-i));
for(k=0;k<=i;k++)
{
ls=m+i-2*k;
tmp=combin(i,i)/(combin(k,k)*combin(i-k,i-k));
work[0*n1+ls]+=tmpd*pow(w0,k)*tmp;
work[1*n1+ls]+=tmpc*pow(w0,k)*tmp;
}
}
for(i=0;i<=n2;i++)
{
d[i]=work[0*n1+i];
c[i]=work[1*n1+i];
}
free(work);
}
}
bilinear(d,c,b,a,n);
}
static double combin(i1,i2)
int i1,i2;
{
int i;
double s;
s=1.0;
if(i2==0) return(s);
for(i=i1;i>(i1-i2);i--)
{s*=i;}
return(s);
}
static void bilinear(d,c,b,a,n)
int n;
double d[],c[],b[],a[];
{
int i,j,n1;
double sum,atmp,scale,*temp;
n1=n+1;
temp=malloc(n1*n1*sizeof(double));
for(j=0;j<=n;j++)
{temp[j*n1+0]=1.0;}
sum=1.0;
for(i=1;i<=n;i++)
{
sum=sum*(double)(n-i+1)/(double)i;
temp[0*n1+i]=sum;
}
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
temp[j*n1+i]=temp[(j-1)*n1+i]-temp[j*n1+i-1]-temp[(j-1)*n1+i-1];
}
for(i=n;i>=0;i--)
{
b[i]=0.0;
atmp=0.0;
for(j=0;j<=n;j++)
{
b[i]=b[i]+temp[j*n1+i]*d[j];
atmp=atmp+temp[j*n1+i]*c[j];
}
scale=atmp;
if(i!=0) a[i]=atmp;
}
for(i=0;i<=n;i++)
{
b[i]=b[i]/scale;
a[i]=a[i]/scale;
}
a[0]=1.0;
free(temp);
}
<< 上一页 [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] ... 下一页 >>