phc15.cpp

基于时域有限差分的遗传算法主程序

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <vector>
#include <iostream>
#include <fstream>
using namespace std;

/* 全局变量 */
double p=0.465e-6;
vector<vector<double> > E;
vector<vector<double> > ER;
double Lstraight=5*p;
double pgrid = 50;
double dx=p/pgrid;
double dz=p/pgrid;

/* 函数定义 */
int round(double,double);
double shape(double,double,double);

double shape(double x,double y,double width)
{
	double z;
    if (x<10)             z=0.5*width*(cos(3.1415*x*0.1)+1)*sqrt(3)*p+0.25*sqrt(3)*p;
	if ((x>9)&&(x<20))    z=0.25*sqrt(3)*p;
    if ((x>19)&&(x<30))   z=0.5*width*(cos((29-x)*3.1415*0.1)+1)*sqrt(3)*p+0.25*sqrt(3)*p;
	return z;
}

int round(double x,double y)
{
	int z;
    z=int (x/y+0.5);
	return z;
}

int txstraight=round(Lstraight,dx);

void main()
{
	int i,j;
    double r=0.3*p;
    double L=40*p;
    double H=7*p;
	double Epsilon=1;
    double EpsilonR=2.8*2.8;
	double dist;
	double distanceL;
	double distanceR;
    double icenter;
    double jcenter;
	double offset;
	double offset1;
	double offset2;
	double width=0.50034;
	double Lcouple=30*p;
	double defectR;
	double defectL;
	double defectcenter=6.25;
	int txcouple=round(Lcouple,dx);
    int tx=round(L,dx);
	int tz=round(H,dz);
	E = vector<vector<double> > (tz,vector<double>(tx,EpsilonR));
    ER = vector<vector<double> > (2*tz,vector<double>(tx,EpsilonR));

	/*char file1[] = "E:\\GA\\SGA\\data\\file1.txt";
	ifstream filedata(file1,ios::in);
	filedata >> width;
	filedata.close();*/

	double Hstraight=(0.5*width+0.5)*sqrt(3)*p;
	int tzstraight=round(Hstraight,dz);
	for (i=0; i<=txstraight-1; i++)
	{
		for(j=0; j<=tz-1; j++)
		{
			if (j>tzstraight-1)
				E[j][i]=Epsilon;
		}
	}
	
    for (i=0; i<=txcouple-1 ; i++)
	{
		for (j=0; j<=tz-1; j++)
		{
			icenter = floor((i+0.5)*dx/p);
		    if(icenter < 10) 
			{
				jcenter=floor((j+0.5)*dz/(0.5*sqrt(3)*p)-offset/(0.5*sqrt(3)*p));
                offset=shape(icenter,jcenter,width);
				offset1=shape(icenter+1,jcenter,width);
			    if((jcenter > -1) && (jcenter < 5))
				{
					if((int)jcenter%2!=0)
						dist=sqrt(pow(((i+0.5)*dx-(icenter+0.5)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
					else
                        dist=10000;
				    if (dist<=r)
						E[j][i+txstraight]=Epsilon;

					if((int)jcenter%2==0)
					{
						distanceL=sqrt(pow(((i+0.5)*dx-icenter*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
						distanceR=sqrt(pow(((i+0.5)*dx-(icenter+1)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset1)),2.0));
					}	
					else
					{
                        distanceL=10000;
					    distanceR=10000;
					}
				    if ((distanceL<=r)||(distanceR<=r))
						E[j][i+txstraight]=Epsilon;
				}
			}
			
            if((icenter > 9) && (icenter < 20))
			{
                jcenter=floor((j+0.5)*dz/(0.5*sqrt(3)*p)-offset/(0.5*sqrt(3)*p));
				offset=shape(icenter,jcenter,width);
				offset1=shape(9.0,jcenter,width);
				offset2=shape(20.0,jcenter,width);
                if((jcenter > -1) && (jcenter < 5))
				{
                    if((int)jcenter%2!=0)
						dist=sqrt(pow(((i+0.5)*dx-(icenter+0.5)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0)); 
                    else
                        dist=10000;
                    if (dist<=r)
					    E[j][i+txstraight]=Epsilon;

				    if((int)jcenter%2==0)
					{
						if(icenter==10) 
						{
							distanceL=sqrt(pow(((i+0.5)*dx-icenter*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset1)),2.0));
							distanceR=sqrt(pow(((i+0.5)*dx-(icenter+1)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0)); 
						}
						else
						{
                            distanceL=10000;
					        distanceR=10000;
						}
				        if ((distanceL<=r)||(distanceR<=r))
						    E[j][i+txstraight]=Epsilon;
						
						if(icenter==19)
						{
							distanceL=sqrt(pow(((i+0.5)*dx-icenter*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
							distanceR=sqrt(pow(((i+0.5)*dx-(icenter+1)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset2)),2.0)); 
						}
						else
						{
                            distanceL=10000;
					        distanceR=10000;
						}
				        if ((distanceL<=r)||(distanceR<=r))
						    E[j][i+txstraight]=Epsilon;
						
						if((icenter>10)&&(icenter<19))
						{
							distanceL=sqrt(pow(((i+0.5)*dx-icenter*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
							distanceR=sqrt(pow(((i+0.5)*dx-(icenter+1)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0)); 
						}
						else
						{
                            distanceL=10000;
					        distanceR=10000;
						}
				        if ((distanceL<=r)||(distanceR<=r))
						    E[j][i+txstraight]=Epsilon;
					}
				}
			}
                    
            if((icenter > 19) && (icenter < 30))
			{
                jcenter=floor((j+0.5)*dz/(0.5*sqrt(3)*p)-offset/(0.5*sqrt(3)*p));
				offset=shape(icenter,jcenter,width);
				offset1=shape(icenter+1,jcenter,width);
                if((jcenter > -1) && (jcenter < 5))
				{
					if((int)jcenter%2!=0)
						dist=sqrt(pow(((i+0.5)*dx-(icenter+0.5)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
					else
                        dist=10000;
                    if (dist<=r)
                        E[j][i+txstraight]=Epsilon;

					if((int)jcenter%2==0)
					{
						distanceL=sqrt(pow(((i+0.5)*dx-icenter*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset)),2.0));
						distanceR=sqrt(pow(((i+0.5)*dx-(icenter+1)*p),2.0)+pow(((j+0.5)*dz-((jcenter+0.5)*0.5*sqrt(3)*p+offset1)),2.0));
					}
					else
					{
                        distanceL=10000;
					    distanceR=10000;
					}
				    if ((distanceL<=r)||(distanceR<=r))
						E[j][i+txstraight]=Epsilon;
				}
            }
		}		
	}
	
	for (i=0; i<=txstraight-1; i++)
	{
		for(j=0; j<=tz-1; j++)
		{
			if (j>tzstraight-1)
				E[j][i+txstraight+txcouple]=Epsilon;
		}
	}

	for (i=0; i<=tx-1; i++)
	{
		for (j=0; j<=tz-1; j++)
		{
			ER[j+tz][i]=E[j][i];
	        ER[tz-j-1][i]=E[j][i];
		}
	}
	
	for(j=0; j<=2*tz-1; j++)
	{
		for(i=0; i<=tx-1; i++)
		{
			defectL=sqrt(pow(((i+0.5)*dx-(defectcenter+5)*p),2.0)+pow(((j+0.5)*dz-7*p),2.0));
            defectR=sqrt(pow(((i+0.5)*dx-(35-defectcenter)*p),2.0)+pow(((j+0.5)*dz-7*p),2.0));
            if (defectL<=r) ER[j][i]=Epsilon;
			if (defectR<=r) ER[j][i]=Epsilon;
		}
			
	}

	char file[100];
    sprintf(file,"E:\\PhC\\PhC data\\file15.txt");
	ofstream eps(file,ios::out);
	for(j=0; j<=2*tz-1; j++)
	{
		for(i=0; i<=tx-1; i++)
			eps<< ER[j][i] << '\t';
		eps<< endl;
	}
	eps.close();
}