?? adifdtd1.cpp
字號:
for (it=0;it<=itmax;it++)
{ //-------------Time interative---------->
if(it%snapStep == 0)
printf("it=%d\n",it);
for(i=0;i<=imax;i++)
{
for(j=0;j<=jmax;j++)
{
if(i==i0 && j==j0 && it<=nstop){/*excitation*/
HZINC1[i][j]=0.5*CFLN*exp(-((it)*tt-t0)*((it)*tt-t0)/(t*t));
HZINC2[i][j]=0.5*CFLN*exp(-((it+0.5)*tt-t0)*((it+0.5)*tt-t0)/(t*t));
}
else
{
HZINC1[i][j]=0.0;
HZINC2[i][j]=0.0;
}
}
}
////////////////////////////////////////////////////////////////////
/*The start of the 1st-updating*/
/*EX field component -- Y_CUT*/
for(i=0;i<=imax;i++)
{
for(j=0;j<=jmax;j++)
{
if(j ==0)
aa[i][j] = -(t2e*t2h)/(yy*yy);
else
aa[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j-1]*yy*yy);
cc[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j]*yy*yy);
bb[i][j] = 1 -aa[i][j] - cc[i][j];
// aa[i][j]=-(t2e*t2h)/(yy*yy);
// bb[i][j]=1.0+2*(t2e*t2h)/(yy*yy);
// cc[i][j]=-(t2e*t2h)/(yy*yy);
}
}
for(i=0;i<imax;i++)
{
for(j=1;j<jmax;j++)
{
r[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*EX[i][j] +
+ (2*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*yy)*(HZ[i][j] - HZ[i][j-1])
- (tt*tt)/(xx*yy*muArray[i][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(EY[i+1][j]-EY[i][j])
+ (tt*tt)/(xx*yy*muArray[i][j-1]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(EY[i+1][j-1]-EY[i][j-1]);
// r[i][j]=EX[i][j]+(t2e/yy)*HZ[i][j]-(t2e/yy)*HZ[i][j-1]
// -((t2e*t2h)/(xx*yy))*(EY[i+1][j]-EY[i][j]) +((t2e*t2h)/(xx*yy))*(EY[i+1][j-1]-EY[i][j-1]);
}
}
for(i=0;i<imax;i++)
{ /*------------------------------------------>*/
for(j=1;j<jmax;j++)
{
if(j==1)
{
bet1 = bb[i][j];
IEX[i][j] = r[i][j]/bet1;
}
else// if (j>=2 && j<=jmax-1)
{
g1[i][j] = cc[i][j-1]/bet1; /*追趕法解3對角*/
bet1 = bb[i][j] - aa[i][j]*g1[i][j];
IEX[i][j] = (r[i][j] - aa[i][j]*IEX[i][j-1])/bet1;
}
}
for(j=jmax-2;j>=1;j--)
{
IEX[i][j]=IEX[i][j]-g1[i][j+1]*IEX[i][j+1];
}
} /*-------------------------------------------->*/
/*EY field component -- Z_CUT*/
for (i=1;i<imax;i++)
{
for (j=0;j<jmax;j++)
{
IEY[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*EY[i][j]
- (2*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*(HZ[i][j]-HZ[i-1][j])/xx;
// IEY[i][j]=EY[i][j]-t2e*(HZ[i][j]-HZ[i-1][j])/xx;
}
}
/*HZ field components*/
for (i=0;i<imax;i++)
{
for (j=0;j<jmax;j++)
{
IHZ[i][j] = HZ[i][j] + HZINC1[i][j]
+ tt/(2*muArray[i][j])*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);
// IHZ[i][j]=HZ[i][j]+HZINC1[i][j] + t2h*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);
}
}
/*end of the 1st-updating*/
///////////////////////////////////////////////////////////////////////
/*start of the 2nd-updating*/
/*EX field component -- Z_CUT*/
for (i=0;i<imax;i++)
{
for (j=1;j<jmax;j++)
{
EX[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*IEX[i][j]
+ (2*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*(IHZ[i][j]-IHZ[i][j-1])/yy;
// EX[i][j]=IEX[i][j]+t2e*(IHZ[i][j]-IHZ[i][j-1])/yy;
}
}
/*EY field component -- X_CUT*/
for(i=0;i<=imax;i++)
{
for(j=0;j<=jmax;j++)
{
if(i==0)
aa[i][j]=-(t2e*t2h)/(xx*xx);
else
aa[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i-1][j]*xx*xx);
cc[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j]*xx*xx);
bb[i][j] = 1 -aa[i][j] - cc[i][j];
// aa[i][j]=-(t2e*t2h)/(xx*xx);
// bb[i][j]=1.0+2*(t2e*t2h)/(xx*xx);
// cc[i][j]=-(t2e*t2h)/(xx*xx);
}
}
for (j=0;j<jmax;j++)
{
for(i=1;i<imax;i++)
{
r[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*IEY[i][j] +
+ (2*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*xx)*(IHZ[i-1][j] - IHZ[i][j])
- (tt*tt)/(xx*yy*muArray[i][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(IEX[i][j+1]-IEX[i][j])
+ (tt*tt)/(xx*yy*muArray[i-1][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(IEX[i-1][j+1]-IEX[i-1][j]);
// r[i][j]=IEY[i][j]-(t2e/xx)*(IHZ[i][j] - IHZ[i-1][j])
// -((t2e*t2h)/(xx*yy))*(IEX[i][j+1]-IEX[i][j]) +((t2e*t2h)/(xx*yy))*(IEX[i-1][j+1]-IEX[i-1][j]);
}
}
for (j=0;j<jmax;j++)//--------------------------------------->
{
for(i=1;i<imax;i++)
{
if(i==1)
{
bet2 = bb[i][j];
EY[i][j] = r[i][j]/bet2;
}
else// if(i>=2&&i<=imax-1)
{
g2[i][j] = cc[i-1][j]/bet2; /*追趕法解3對角*/
bet2 = bb[i][j] - aa[i][j]*g2[i][j];
EY[i][j] = (r[i][j] - aa[i][j]*EY[i-1][j])/bet2;
}
}
for(i=imax-2;i>=1;i--)
{
EY[i][j]=EY[i][j]-g2[i+1][j]*EY[i+1][j];
}
} /*-------------------------------------------->*/
/*HZ field components*/
for (i=0;i<imax;i++)
{
for (j=0;j<jmax;j++)
{
HZ[i][j] = IHZ[i][j] + HZINC2[i][j]
+ tt/(2*muArray[i][j])*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);
// HZ[i][j]=IHZ[i][j]+HZINC2[i][j]
// +t2h*((IEX[i][j+1]-IEX[i][j])/yy
// -(EY[i+1][j]-EY[i][j])/xx);
}
}
/*end of the 2nd-updating*/
/*
char *strEXtemp = "EX_";
char *strEYtemp = "EY_";
char *strHZtemp = "HZ_";
char stimulusFilename[30];
char addString[40];
if(fmod(it, snapStep)< 1e-8)
{
strcpy(stimulusFilename, strEXtemp);
sprintf(addString, "%d.dat", it);
strcat(stimulusFilename, addString);
fp3 = fopen(stimulusFilename, "w");
for(int ik = 0; ik < imax; ik++)
{
for(int jk =0; jk < jmax; jk++)
{
fprintf(fp3,"%15.7f ", EX[ik][jk]);
}
fprintf(fp3,"\n");
}
fclose(fp3);
strcpy(stimulusFilename, strEYtemp);
sprintf(addString, "%d.dat", it);
strcat(stimulusFilename, addString);
fp3 = fopen(stimulusFilename, "w");
for( ik = 0; ik < imax; ik++)
{
for(int jk =0; jk < jmax; jk++)
{
fprintf(fp3,"%15.7f ", EY[ik][jk]);
}
fprintf(fp3,"\n");
}
fclose(fp3);
strcpy(stimulusFilename, strHZtemp);
sprintf(addString, "%d.dat", it);
strcat(stimulusFilename, addString);
fp3 = fopen(stimulusFilename, "w");
for( ik = 0; ik < imax; ik++)
{
for(int jk =0; jk < jmax; jk++)
{
fprintf(fp3,"%15.7f ", HZ[ik][jk]);
}
fprintf(fp3,"\n");
}
fclose(fp3);
}
*/
fprintf(fp1,"%15.6f, %15.7f, %15.7f, %15.7f, %15.7f, %15.7f\n",it*tt*(10e9),EX[nobsX][nobsY],EY[nobsX][nobsY],
HZ[nobsX][nobsY], HZINC1[i0][j0], HZINC2[i0][j0]);
/*recording the field component*/
H1[it]=HZ[int(imax*0.9)][int(jmax*0.5)];
} //-------------Time interative---------->
printf("Time used before FT = %5d (seconds)\n",(time(NULL)-myt));
for(k=1;k<1499;k++){
xre1=xim1=0.0;
ff[k]=k/(tt*itmax);
w=2.0*pi*ff[k];
for(i=0;i<=itmax;i++)
{
xre1=xre1+H1[i]*cos(w*tt*i);
xim1=xim1-H1[i]*sin(w*tt*i);
}
rr1[k]=(xre1*xre1+xim1*xim1);
fprintf(fp2,"%15.7f, %15.7f\n",ff[k]/10e+9, rr1[k]);
}
fclose (fp1);
fclose (fp2);
/*stuff related to C++*/
// delete EX, EY, HZ
for(i=0;i<=imax;i++){
delete [] EX[i];
}
delete [] EX;
for(i=0;i<=imax;i++){
delete [] EY[i];
}
delete [] EY;
for(i=0;i<=imax;i++){
delete [] HZ[i];
}
delete [] HZ;
for(i=0;i<=imax;i++){
delete [] IEX[i];
}
delete [] IEX;
for(i=0;i<=imax;i++){
delete [] IEY[i];
}
delete [] IEY;
for(i=0;i<=imax;i++){
delete [] IHZ[i];
}
delete [] IHZ;
// delete HZINC1, HZINC2
for(i=0;i<=imax;i++)
delete [] HZINC1[i];
delete [] HZINC1;
for(i=0;i<=imax;i++)
delete [] HZINC2[i];
delete [] HZINC2;
/////////////////////
for(i=0;i<=imax;i++){
delete [] thetaArray[i];
}
delete [] thetaArray;
for(i=0;i<=imax;i++){
delete [] muArray[i];
}
delete [] muArray;
for(i=0;i<=imax;i++){
delete [] epsArray[i];
}
delete [] epsArray;
/////////////////////
for(i=0;i<=imax;i++){
delete [] aa[i];
}
delete [] aa;
for(i=0;i<=imax;i++){
delete [] bb[i];
}
delete [] bb;
for(i=0;i<=imax;i++){
delete [] cc[i];
}
delete [] cc;
for(i=0;i<=imax;i++){
delete [] r[i];
}
delete [] r;
for(i=0;i<=imax;i++){
delete [] g1[i];
}
delete [] g1;
for(i=0;i<=imax;i++){
delete [] g2[i];
}
delete [] g2;
delete [] H1;
printf("Total CPU Time used = %5d (seconds)\n",(time(NULL)-myt));
return 0;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -