nl=m0;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx0;
         r=r+tt[l][ml]*xez2[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempz[i1]=tempz[i1]+temp;
       nl++;} }
   for(i=xs;i<=xf;i++)
     { if(zz==Z2) ez2[i][yy][zz]=tempz[i+30];
       else ez1[i][yy][zz]=tempz[i]; }
   for (i=xs;i<=xf;i++)
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempx1[i1]=0;
     nl=m;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx;
         r=r+tt[l][ml]*xex1[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempx1[i1]=tempx1[i1]+temp;
       nl++;} }
   for (i=xs;i<=xf;i++)
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempx2[i1]=0;
     nl=m0;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx0;
         r=r+tt[l][ml]*xex1[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempx2[i1]=tempx2[i1]+temp;
       nl++;} }
   for(i=xs;i<=xf;i++)
     { if(zz==Z2) ex2[i][yy][zz]=(tempx1[i+30]+tempx2[i+30])/2;
       else ex1[i][yy][zz]=(tempx1[i]+tempx2[i])/2; }
   for(i=0;i<=X1;i++)
    for(j=0;j<=14;j++)
     for(nl=6;nl>=1;nl--)
     { i1=j+tx;
       nl0=nl+m;
       xey1[i][i1][nl0]=xey1[i][i1][nl0-1];
       xex1[i][i1][nl0]=xex1[i][i1][nl0-1]; }
   for(i=0;i<=X1;i++)
    for(j=0;j<=14;j++)
     for(nl=6;nl>=1;nl--)
     { i1=j+tx0;
       nl0=nl+m0;
       xez2[i][i1][nl0]=xez2[i][i1][nl0-1];
       xex2[i][i1][nl0]=xex2[i][i1][nl0-1]; }
   for(i=xs;i<=xf;i++)
    for(j=0;j<=14;j++)
    { i1=j+tx;
      if(yy==0) i2=yy+i;
      else i2=yy-i;
      if (zz==Z2){xey1[i+30][i1][m]=ey2[i][i2][zz];
                  xex1[i+30][i1][m]=ex2[i][i2][zz];}
      else {xey1[i][i1][m]=ey1[i][i2][zz];
            xex1[i][i1][m]=ex1[i][i2][zz];}    }
   for(i=xs;i<=xf;i++)
    for(j=0;j<=14;j++)
    { i1=j+tx0;
      if(yy==0) i2=yy+i;
      else i2=yy-i;
      if (zz==Z2){xez2[i+30][i1][m0]=ez2[i][i2][zz];
                  xex2[i+30][i1][m0]=ex2[i][i2][zz];}
      else {xez2[i][i1][m0]=ez1[i][i2][zz];
            xex2[i][i1][m0]=ex1[i][i2][zz];}    }

 }

///////////////////////////////////////////////////////////////////////////////
 void  fourier(float e[X2+1][Y+1][rn+1])
  {  int    n,k,i,j;
     float  gk,gre[X2+1][Y+1][4],gim[X2+1][Y+1][4],a[rn+1];
     gk=8.5/df;
     for (k=0;k<=3;k++)
      for(i=0;i<=X2;i++)
       for(j=0;j<=Y;j++)
          {gre[i][j][k]=0;
           gim[i][j][k]=0;}
     for(k=0;k<=3;k++)
      { for(n=0;n<=rn;n++)
        { a[n]=(2*pi*gk*n)/rn; }
        for(i=0;i<=X2;i++)
         for(j=0;j<=Y;j++)
          for(n=0;n<=rn;n++)
            {gre[i][j][k]=gre[i][j][k]+e[i][j][n]*cos(a[n]);
             gim[i][j][k]=gim[i][j][k]+e[i][j][n]*sin(a[n]); }
        gk=gk+1/df; }
     if (e==rex)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fex[i][j][0]=gre[i][j][1];
                  fex[i][j][1]=gim[i][j][1];}    }
     if (e==rey)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fey[i][j][0]=gre[i][j][1];
                  fey[i][j][1]=gim[i][j][1];}    }
     if (e==rhx)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fhx[i][j][0]=gre[i][j][1];
                  fhx[i][j][1]=gim[i][j][1];}    }
     if (e==rhy)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fhy[i][j][0]=gre[i][j][1];
                  fhy[i][j][1]=gim[i][j][1];}    }
  }

//////////////////////////////////////////////////////////////////////////////////////////
 void sss()
{ int a,b,c0;
  float k,gr1,gi1,gr2,gi2,temp;
  float gr[X2+1][Y+1][46],gi[X2+1][Y+1][46];
  k=2*pi/ld;
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fey[a][b][0]-ita*fhx[a][b][0]*cos(c0*dc);
       gi1=fey[a][b][1]-ita*fhx[a][b][1]*cos(c0*dc);
       temp=k*((-1)*(b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er1[a][b][c0+45]=gi[a][b][c0+45]*(k/(4*pi));
       Ei1[a][b][c0+45]=(-1)*gr[a][b][c0+45]*(k/(4*pi));  }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fey[a][b][0]+ita*fhx[a][b][0]*cos(c0*dc);
       gi1=-fey[a][b][1]+ita*fhx[a][b][1]*cos(c0*dc);
       temp=k*((b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er1[a][b][c0+45]=gi[a][b][c0]*(k/(4*pi));
       Ei1[a][b][c0+45]=(-1)*gr[a][b][c0]*(k/(4*pi));}
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se1[c0][0]=se1[c0][0]+Er1[a][b][c0];
       se1[c0][1]=se1[c0][1]+Ei1[a][b][c0];}
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fex[a][b][0]*cos(c0*dc)+ita*fhy[a][b][0];
       gi1=fex[a][b][1]*cos(c0*dc)+ita*fhy[a][b][1];
       temp=k*((-1)*(b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er2[a][b][c0+45]=(-1)*gi[a][b][c0+45]*(k/(4*pi));
       Ei2[a][b][c0+45]=gr[a][b][c0+45]*(k/(4*pi));    }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fex[a][b][0]*cos(c0*dc)-ita*fhy[a][b][0];
       gi1=-fex[a][b][1]*cos(c0*dc)-ita*fhy[a][b][1];
       temp=k*((b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er2[a][b][c0+45]=(-1)*gi[a][b][c0]*(k/(4*pi));
       Ei2[a][b][c0+45]=gr[a][b][c0]*(k/(4*pi));       }
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se2[c0][0]=se2[c0][0]+Er2[a][b][c0];
       se2[c0][1]=se2[c0][1]+Ei2[a][b][c0];}
  for(c0=0;c0<=90;c0++)
   {  Me1[c0]=sqrt(pow(se1[c0][0],2)+pow(se1[c0][1],2));
      Me2[c0]=sqrt(pow(se2[c0][0],2)+pow(se2[c0][1],2)); }
}

///////////////////////////////////////////////////////////////////////////////////////////
 void sss1()
{ int a,b,c0;
  float k,gr1,gi1,gr2,gi2,temp;
  float gr[X2+1][Y+1][46],gi[X2+1][Y+1][46];
  k=2*pi/ld;
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fex[a][b][0]-ita*fhy[a][b][0]*cos(c0*dc);
       gi1=-fex[a][b][1]-ita*fhy[a][b][1]*cos(c0*dc);
       temp=k*((-1)*(a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er3[a][b][c0+45]=gi[a][b][c0+45]*(k/(4*pi));
       Ei3[a][b][c0+45]=(-1)*gr[a][b][c0+45]*(k/(4*pi));  }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fex[a][b][0]+ita*fhy[a][b][0]*cos(c0*dc);
       gi1=fex[a][b][1]+ita*fhy[a][b][1]*cos(c0*dc);
       temp=k*((a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er3[a][b][c0+45]=gi[a][b][c0]*(k/(4*pi));
       Ei3[a][b][c0+45]=(-1)*gr[a][b][c0]*(k/(4*pi));}
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se3[c0][0]=se3[c0][0]+Er3[a][b][c0];
       se3[c0][1]=se3[c0][1]+Ei3[a][b][c0];}
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fey[a][b][0]*cos(c0*dc)-ita*fhx[a][b][0];
       gi1=fey[a][b][1]*cos(c0*dc)-ita*fhx[a][b][1];
       temp=k*((-1)*(a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er4[a][b][c0+45]=(-1)*gi[a][b][c0+45]*(k/(4*pi));
       Ei4[a][b][c0+45]=gr[a][b][c0+45]*(k/(4*pi));    }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fey[a][b][0]*cos(c0*dc)+ita*fhx[a][b][0];
       gi1=-fey[a][b][1]*cos(c0*dc)+ita*fhx[a][b][1];
       temp=k*((a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er4[a][b][c0+45]=(-1)*gi[a][b][c0]*(k/(4*pi));
       Ei4[a][b][c0+45]=gr[a][b][c0]*(k/(4*pi));       }
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se4[c0][0]=se4[c0][0]+Er4[a][b][c0];
       se4[c0][1]=se4[c0][1]+Ei4[a][b][c0];}
  for(c0=0;c0<=90;c0++)
   {  Me3[c0]=sqrt(pow(se3[c0][0],2)+pow(se3[c0][1],2));
      Me4[c0]=sqrt(pow(se4[c0][0],2)+pow(se4[c0][1],2)); }
}

//////////////////////////////////////////////////////////////////////////////////////////////////
void fourier1(float e[rn+1])
{ int n,k;
  float gk,gr[4],gi[4],a[rn+1];
  gk=8.5/df;
  for (k=0;k<=3;k++)
   { for(n=0;n<=rn;n++)
      { a[n]=(2*pi*gk*n)/rn;  }
     gr[k]=0;gi[k]=0;
     for(n=0;n<=rn;n++)
      { gr[k]=gr[k]+e[n]*cos(a[n]);
        gi[k]=gi[k]+e[n]*sin(a[n]);}
     gk=gk+1/df; }
  if (e==Vs) {Vsw[0]=gr[1];
              Vsw[1]=gi[1];}
  if (e==Is) {Isw[0]=gr[1];
              Isw[1]=gi[1];}
}

//////////////////////////////////////////////////////////////////////////////////////////////
void mo(float e[X2+1][Y+1][2])
{FILE *fp;
 int i,j;
 float mod[X2+1][Y+1];
 for(i=0;i<=X2;i++)
  for(j=0;j<=Y;j++)
  { mod[i][j]=sqrt(pow(e[i][j][0],2)+pow(e[i][j][1],2)); }
 if (e==fex){ fp=fopen("mex.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mex[i][j]=mod[i][j];
                 fprintf(fp,"mex(%d,%d)=%.12f;\n",i+1,j+1,mex[i][j]);}
             fprintf(fp,"mesh(mex);");
             fclose(fp);}
 if (e==fey){ fp=fopen("mey.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mey[i][j]=mod[i][j];
                 fprintf(fp,"mey(%d,%d)=%.12f;\n",i+1,j+1,mey[i][j]);}
             fprintf(fp,"mesh(mey);");
             fclose(fp);}
 if (e==fhx){ fp=fopen("mhx.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mhx[i][j]=mod[i][j];
                 fprintf(fp,"mhx(%d,%d)=%.12f;\n",i+1,j+1,mhx[i][j]);}
             fprintf(fp,"mesh(mhx);");
             fclose(fp);}
 if (e==fhy){ fp=fopen("mhy.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mhy[i][j]=mod[i][j];
                 fprintf(fp,"mhy(%d,%d)=%.12f;\n",i+1,j+1,mhy[i][j]);}
             fprintf(fp,"mesh(mhy);");
             fclose(fp);}
}

//////////////////////////////////////////////////////////////////////////////
//輸出數據
void putdata() {
  FILE  *fp14,*fp15,*fp16,*fp17,*fp18;
 int n;
  float tmp;
 fp14=fopen("pin.txt","w");
 pin[0]=Vsw[0]*Isw[0]+Vsw[1]*Isw[1];
 pin[1]=-Vsw[0]*Isw[1]+Vsw[1]*Isw[0];
 Mp=sqrt(pow(pin[0],2)+pow(pin[1],2));
 fprintf(fp14,"pin=%.12f+j%.12f\nMode=%.12f\n",pin[0],pin[1],Mp);
 fclose(fp14);
 for(n=0;n<=90;n++)
  { tmp=fabs(pow(Me1[n],2)*pow(dr,2)/(30*Mp));
    G1[n]=10*log10(tmp); }
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me2[n],2)*pow(dr,2)/(30*Mp));
   G2[n]=10*log10(tmp); }
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me3[n],2)*pow(dr,2)*cos((n0-45)*dc)/(30*Mp));
   G3[n]=10*log10(tmp);}
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me4[n],2)*pow(dr,2)/(30*Mp));
   G4[n]=10*log10(tmp);}
 fp15=fopen("G1.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp15,"y(%d)=%.12f;\n",n+1,G1[n]);
 fprintf(fp15,"plot(y);");
 fclose(fp15);
 fp16=fopen("G2.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp16,"y(%d)=%.12f;\n",n+1,G2[n]);
 fprintf(fp16,"plot(y);");
 fclose(fp16);
 fp17=fopen("G3.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp17,"y(%d)=%.12f;\n",n+1,G3[n]);
 fprintf(fp17,"plot(y);");
 fclose(fp17);
 fp18=fopen("G4.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp18,"y(%d)=%.12f;\n",n+1,G4[n]);
 fprintf(fp18,"plot(y);");
 fclose(fp18);
}