simulateview.cpp

运用扫描法求设计一个最佳比例的聚光腔

//			AfxBeginThread(TestThread, (LPVOID)param);	
			m_bStop =false;
			DWORD dwThreadID;
			DWORD (WINAPI CSimulateView:: *p1)() =SimulateProc;
			DWORD (WINAPI *p2)(LPVOID) = NULL;		
			memcpy(&p2, &p1, sizeof(p1));
			// Be sure to pass "this" parameter to the thread procedure.
			m_hThread = ::CreateThread(NULL, 0, p2, this, 0, &dwThreadID);
			if(m_hThread == NULL)
			{
				MessageBox("创建线程时出现错误！","嗅探",MB_ICONINFORMATION);
				return;	
			}
		}	

	//判断所有线程是否结束，否则空循环等待
	while(m_nCurThreadNum >0)
	{};*/
}


//对于当前的长半轴和短半轴随机生成已经给出的光子数并模拟
DWORD WINAPI CSimulateView::SingleSimulateProc()
{
	Point pCutPoint;   //交点
	Point pTempPoint;  //临时交点
	Line pReflectLine; //反射直线
	int nReflectCount; //已反射次数
	bool IsDie=false;  //光子是否消灭
	bool IsSorb=false; //光子是否被吸收
	
	//清除最大记录
	m_rMaxRecord.fGoodRate=0;

	m_Rand.SetTimerSeed();

	for(m_fCurLAxis=m_fLAxisFrom;m_fCurLAxis<=80;m_fCurLAxis+=m_fLAxisAdd)
		for(m_fCurSAxis=Max(m_fSAxisFrom,(float)sqrt(2*m_fCurLAxis*m_fObjRadius-m_fObjRadius*m_fObjRadius));m_fCurSAxis<m_fCurLAxis;m_fCurSAxis+=m_fSAxisAdd)
		{
			m_nGoodPhotonCount=0;
			m_nBadPhotonCount=0;
			m_nDeadPhotonCount=0;

			m_fFocus=this->GetFocus(this->m_fCurLAxis);
			
			for(int i=1;i<=this->m_nAllPhoton;i++)
			{
				IsDie=false;
				
				if(m_nObjSort==0)//球面
				{
					//-----------产生随机光子并发射之--------------					
					pReflectLine.dPointDirct.dx=(float)m_Rand.DrawRandomNumber();
					pReflectLine.dPointDirct.dy=(float)m_Rand.DrawRandomNumber();
					pReflectLine.dPointDirct.dz=(float)m_Rand.DrawRandomNumber();
					
					pReflectLine.pPassPoint.x=this->GetSAxis(this->m_fCurLAxis);
					pReflectLine.pPassPoint.y=(float)m_Rand.DrawRandomNumber()/100*2*m_fCurSAxis*m_fCurSAxis/m_fCurLAxis-m_fCurSAxis*m_fCurSAxis/m_fCurLAxis;
					pReflectLine.pPassPoint.z=(float)m_Rand.DrawRandomNumber()/100*2*m_fCurSAxis*m_fCurSAxis/m_fCurLAxis-m_fCurSAxis*m_fCurSAxis/m_fCurLAxis;
					//---------------------------------------------
				}
				else
				{
					//-----------产生随机光子并发射之--------------					
					pReflectLine.dPointDirct.dx=(float)m_Rand.DrawRandomNumber();
					pReflectLine.dPointDirct.dy=(float)m_Rand.DrawRandomNumber();
					pReflectLine.dPointDirct.dz=(float)m_Rand.DrawRandomNumber();
			
					pReflectLine.pPassPoint.x=(float)m_Rand.DrawRandomNumber()/100*2*this->m_fObjRadius+this->m_fFocus-this->m_fObjRadius;
					pReflectLine.pPassPoint.y=(float)m_Rand.DrawRandomNumber()/100*2*this->m_fObjRadius-this->m_fObjRadius;
					pReflectLine.pPassPoint.z=(float)m_Rand.DrawRandomNumber()/100*2*this->m_fObjRadius-this->m_fObjRadius;
					//---------------------------------------------
				}
					
				//-----------------求出交点--------------------
				pCutPoint=this->GetCutPoint(pReflectLine,-1);
				//---------------------------------------------
				
				pReflectLine.pPassPoint.x=pCutPoint.x;
				pReflectLine.pPassPoint.y=pCutPoint.y;
				pReflectLine.pPassPoint.z=pCutPoint.z;

				//光子在没有被光电管接收或被物体消灭之前不停的在椭球内反射
				for(nReflectCount=0;nReflectCount<=this->m_nReflectCount;nReflectCount++)
				{				
					//得出反射直线，反射直线的方程：(x-x1)/x2=(y-y1)/y2=(z-z1)/z2
					pReflectLine=this->GetReflectLine(pReflectLine);
					pTempPoint=pCutPoint;

					//求出反射直线与椭球的交点
					pCutPoint=this->GetCutPoint(pReflectLine,1);
					if(this->IsEven(pCutPoint,pTempPoint))
						pCutPoint=this->GetCutPoint(pReflectLine,-1);
					
					switch(this->IsCutWithPipe(pReflectLine))
					{
					case -1://不与光电管相交
						break;
					case 0://光子消灭
						IsDie=true;						
						break;
					case 1://光子被光电管接收
						IsSorb=true;
						m_nGoodPhotonCount+=1;
						break;
					default:
						break;
					}
					
					//是否与物体相交
					if(this->IsCutWithObj(pReflectLine,m_nObjSort))
						IsDie=true;
					
					if(IsDie)
						m_nBadPhotonCount+=1;

					if(IsDie || IsSorb)
						break;									
				}
			}

			//--------------------------插入记录-----------------------------------
			m_nDeadPhotonCount=m_nAllPhoton-m_nGoodPhotonCount-m_nBadPhotonCount;

			m_fGoodRate=(float)this->m_nGoodPhotonCount/this->m_nAllPhoton;
			m_fBadRate=(float)this->m_nBadPhotonCount/this->m_nAllPhoton;
			m_fDeadRate=(float)this->m_nDeadPhotonCount/this->m_nAllPhoton;

			m_Record.fLAxis=this->m_fCurLAxis;
			m_Record.fSAxis=this->m_fCurSAxis;
			m_Record.nBadPhotonCount=this->m_nBadPhotonCount;
			m_Record.nGoodPhotonCount=this->m_nGoodPhotonCount;
			m_Record.nDeadPhotonCount=this->m_nDeadPhotonCount;
			m_Record.fGoodRate=this->m_fGoodRate;
			m_Record.fBadRate=this->m_fBadRate;
			m_Record.fDeadRate=this->m_fDeadRate;

//			InsertRecord(m_Record);			
			//---------------------------------------------------------------------


			//--------------------------记录最佳效果-------------------------------
			if(m_fGoodRate>m_rMaxRecord.fGoodRate)
			{
				m_rMaxRecord.fLAxis=this->m_fCurLAxis;
				m_rMaxRecord.fSAxis=this->m_fCurSAxis;
				m_rMaxRecord.nIndex=m_pList->GetItemCount();
				m_rMaxRecord.nGoodPhotonCount=this->m_nGoodPhotonCount;
				m_rMaxRecord.fGoodRate=this->m_fGoodRate;
				m_rMaxRecord.fBadRate=this->m_fBadRate;
				m_rMaxRecord.fDeadRate=this->m_fDeadRate;
				m_rMaxRecord.nBadPhotonCount=this->m_nBadPhotonCount;
				m_rMaxRecord.nDeadPhotonCount=this->m_nDeadPhotonCount;
				m_rMaxRecord.nGoodPhotonCount=this->m_nGoodPhotonCount;
//				m_pSimulateViewList->rMaxRecord=m_rMaxRecord;
			}
			//--------------------------------------------------------------------
		}

	return 0;
}


float CSimulateView::Min(float a, float b)
{
	if(a<=b)
		return a;
	else
		return b;
}

float CSimulateView::Max(float a, float b)
{
	if(a>=b)
		return a;
	else
		return b;
}


Point CSimulateView::GetCutPoint(Line l0,int nID)
{
	Point pPoint;
	float a,b,c;//椭球的长短半轴(50<=a<=80):椭圆方程：x*x/a^2+y*y/b^2+z*z/b^2=1
	float px0,py0,pz0,dx0,dy0,dz0; //已知直线：(x-px0)/dx0=(y-py0)/dy0=(z-pz0)/dz0=k
	float cx,cy,cz;//已知直线l与椭球的交点

	a=this->m_fCurLAxis;
	b=this->m_fCurSAxis;
	c=this->m_fFocus;

	px0=l0.pPassPoint.x;
	py0=l0.pPassPoint.y;
	pz0=l0.pPassPoint.z;
	dx0=l0.dPointDirct.dx;
	dy0=l0.dPointDirct.dy;
	dz0=l0.dPointDirct.dz;

	//椭球与已知直线联立方程得：(a2*a2/a^2+b2^2/b^2+c2^2/b^2)k1^2+2*(a2*a1/a^2+b2*b1/b^2+c2*c1/b^2)+(a1^2/a^2+b1^2/b^2+c1^2/b^2-1)=0
	float m,n,p,k;
	m=dx0*dx0/(a*a)+dy0*dy0/(b*b)+dz0*dz0/(b*b);
	n=dx0*px0/(a*a)+dy0*py0/(b*b)+dz0*pz0/(b*b);
	p=px0*px0/(a*a)+py0*py0/(b*b)+pz0*pz0/(b*b)-1;

	if(nID==1)
		k=(float)(-n+sqrt(n*n-m*p))/m;
	else
		if(nID==-1)
			k=(float)(-n-sqrt(n*n-m*p))/m;

	//求出交点
	cx=dx0*k+px0;
	cy=dy0*k+py0;
	cz=dz0*k+pz0;

	pPoint.x=cx;
	pPoint.y=cy;
	pPoint.z=cz;
	
	return pPoint;
}

Line CSimulateView::GetReflectLine(Line l0)
{
	float a,b,c;//椭球的长短半轴(50<=b<=80):椭圆方程：x*x/a^2+y*y/b^2+z*z/b^2=1
	float cx0,cy0,cz0,dx0,dy0,dz0; //已知直线：(x-cx0)/dx0=(y-cy0)/dy0=(z-cz0)/dz0=k
	float dx1,dy1,dz1;//对称轴：a*a*(x-cx0)/cx0=b*b*(y-cy0)/cy0=b*b*(z-cz0)/cz0

	Line lRet;     //返回的直线

	a=this->m_fCurLAxis;
	b=this->m_fCurSAxis;
	c=this->m_fFocus;

	cx0=l0.pPassPoint.x;
	cy0=l0.pPassPoint.y;
	cz0=l0.pPassPoint.z;
	dx0=l0.dPointDirct.dx;
	dy0=l0.dPointDirct.dy;
	dz0=l0.dPointDirct.dz;

	dx1=cx0/(a*a);
	dy1=cy0/(b*b);
	dz1=cz0/(b*b);

	float k;

	k=2*(dx1*dx0+dy1*dy0+dz1*dz0)/(dx1*dx1+dy1*dy1+dz1*dz1);
	lRet.dPointDirct.dx=dx1*k-dx0;
	lRet.dPointDirct.dy=dy1*k-dy0;
	lRet.dPointDirct.dz=dz1*k-dz0;
	lRet.pPassPoint.x=cx0;
	lRet.pPassPoint.y=cy0;
	lRet.pPassPoint.z=cz0;

	return lRet;
}

float CSimulateView::GetSAxis(float fLAxis)
{
	return (float)sqrt(1.41*this->m_fCurLAxis);
}

float CSimulateView::GetFocus(float fLAxis)
{
	return (float)sqrt(this->m_fCurLAxis*this->m_fCurLAxis-this->m_fCurSAxis*this->m_fCurSAxis);
}


bool CSimulateView::IsEven(Point p1, Point p2)
{
	if(p1.x==p2.x&&p1.y==p2.y&&p1.z==p2.z)
		return true;
	else
		return false;
}


bool CSimulateView::IsCutWithObj(Line l0,int nObjSort)
{
	float cx0=0,cy0=0,cz0=0,dx0=0,dy0=0,dz0=0; //已知直线：(x-cx0)/dx0=(y-cy0)/dy0=(z-cz0)/dz0=k
	float r,c;//物体：(x-c)^2+y^2+z^2=r^2;
	float m,n,p;

	r=this->m_fObjRadius;
	c=this->m_fFocus;

	if(nObjSort==0)//球面
	{
		m=dy0*dy0+dz0*dz0;
		n=dy0*cy0+dz0*cz0;
		p=cy0*cy0+cz0*cz0-r*r;
		if(m*n-p*p<0)
			return false;
		else
			return true;
	}

	cx0=l0.pPassPoint.x;
	cy0=l0.pPassPoint.y;
	cz0=l0.pPassPoint.z;
	dx0=l0.dPointDirct.dx;
	dy0=l0.dPointDirct.dy;
	dz0=l0.dPointDirct.dz;


	m=dx0*dx0+dy0*dy0+dz0*dz0;
	n=(cx0-c)*dx0+cy0*dy0+cz0*dz0;
	p=(cx0-c)*(cx0-c)+cy0*cy0+cz0*cz0-r*r;

	if(n*n-m*p<0)
		return false;
	else
		return true;
}

//返回-1表示不与光电管相交，返回0表示光子消灭，返回1表示光子被光电管接收
int CSimulateView::IsCutWithPipe(Line l0)
{
	float cx0,cy0,cz0,dx0,dy0,dz0; //已知直线：(x-cx0)/dx0=(y-cy0)/dy0=(z-cz0)/dz0=k
	float a,c;//光电倍增管方程：y*y+z*z=1.41*1.41 -a<=x<=-c
	float cx1,cx2;

	a=this->m_fCurLAxis;
	c=this->m_fFocus;

	cx0=l0.pPassPoint.x;
	cy0=l0.pPassPoint.y;
	cz0=l0.pPassPoint.z;
	dx0=l0.dPointDirct.dx;
	dy0=l0.dPointDirct.dy;
	dz0=l0.dPointDirct.dz;

	float m,n,p,k1,k2;

	m=dy0*dy0+dz0*dz0;
	n=dy0*cy0+dz0*cz0;
	p=(float)(cy0*cy0+cz0*cz0-1.41*1.41);

	if(n*n-m*p<0)
		return -1;


	//-------------与形式管有交点----------------
	//形式管是光电管的延伸
	k1=(float)(-n+sqrt(n*n-m*p))/m;
	k2=(float)(-n-sqrt(n*n-m*p))/m;
	cx1=dx0*k1+cx0;
	cx2=dx0*k2+cx0; 
	if(cx0<-c)//此时光子不可能被光电管接收
	{
		if(cx1>-c && cx2>-c)
			return -1;
		else
			return 0;
	}
	else
	{
		if(cx1>-c && cx2>-c)
			return -1;
		else
		{
	///		if(cx1<-c && cx2<-c)
	//			return 0;
	//		else
	//			return 1;
			if((cx1+c)*(cx2+c)<=0)
				return 1;
			else
				return 0;
		}
	}
	//-----------------------------------------
}

void CSimulateView::OnDblclk(NMHDR* pNMHDR, LRESULT* pResult) 
{
	// TODO: Add your control notification handler code here
	this->OnProperty();

	*pResult = 0;
}

void CSimulateView::OnProperty() 
{
	// TODO: Add your command handler code here
	int nIndex;
	nIndex=m_pList->GetItemCount();
	if(nIndex==0)
		return;

	CString strNum,strLAxis,strSAxis,strGoodPhoton,strGoodRate,strBadPhoton,strBadRate,strDeadPhoton,strDeadRate;

	strNum=m_pList->GetItemText((int)m_nCurSel - 1,0);
	strLAxis=m_pList->GetItemText((int)m_nCurSel - 1,1);
	strSAxis=m_pList->GetItemText((int)m_nCurSel - 1,2);
	strGoodPhoton=m_pList->GetItemText((int)m_nCurSel - 1,3);
	strGoodRate=m_pList->GetItemText((int)m_nCurSel - 1,4);
	strBadPhoton=m_pList->GetItemText((int)m_nCurSel - 1,5);
	strBadRate=m_pList->GetItemText((int)m_nCurSel - 1,6);
	strDeadPhoton=m_pList->GetItemText((int)m_nCurSel - 1,7);
	strDeadRate=m_pList->GetItemText((int)m_nCurSel - 1,8);

	CPropertyDlg dlg;
	dlg.DispData(strNum,strLAxis,strSAxis,strGoodPhoton,strGoodRate,strBadPhoton,strBadRate,strDeadPhoton,strDeadRate);
	dlg.DoModal();
}

void CSimulateView::AnalyseData(float fError,CSimulateView* p)
{
	float fMaxGoodRate=m_rMaxRecord.fGoodRate;
	float fGoodRate;
	CString strGoodRate;
	int i=0,nIndex;
	char* pGoodRate;
	RecordString rs;
	int nCount=m_pList->GetItemCount();

	for(i=0;i<nCount;i++)
	{
		strGoodRate=m_pList->GetItemText(i,4);
		pGoodRate=(LPSTR)(LPCTSTR)strGoodRate.Left(strGoodRate.GetLength()-1);
		fGoodRate=atoi(pGoodRate);
		if(abs(fGoodRate-fMaxGoodRate)<=fError)
		{
			nIndex=p->m_pList->GetItemCount();
			rs=this->CopyRecord(nIndex);
			p->InsertRecord(rs);
		}
	}
}

RecordString CSimulateView::CopyRecord(int nIndex)
{
	RecordString r;

	r.strLAxis=m_pList->GetItemText(nIndex,0);
	r.strSAxis=m_pList->GetItemText(nIndex,1);
	r.strGoodPhotonCount=m_pList->GetItemText(nIndex,2);
	r.strBadPhotonCount=m_pList->GetItemText(nIndex,3);
	r.strDeadPhotonCount=m_pList->GetItemText(nIndex,4);
	r.strGoodRate=m_pList->GetItemText(nIndex,5);
	r.strBadRate=m_pList->GetItemText(nIndex,6);
	r.strDeadRate=m_pList->GetItemText(nIndex,7);

	return r;
}

void CSimulateView::InsertRecord(RecordString rs)
{
	this->m_pList->InsertItem(m_nIndex,rs.strIndex);
	this->m_pList->SetItemText(m_nIndex,1,rs.strLAxis);
	this->m_pList->SetItemText(m_nIndex,2,rs.strSAxis);
	this->m_pList->SetItemText(m_nIndex,3,rs.strGoodPhotonCount);
	this->m_pList->SetItemText(m_nIndex,4,rs.strGoodRate);

	this->m_pList->SetItemText(m_nIndex,5,rs.strBadPhotonCount);
	this->m_pList->SetItemText(m_nIndex,6,rs.strBadRate);

	this->m_pList->SetItemText(m_nIndex,7,rs.strDeadPhotonCount);
	this->m_pList->SetItemText(m_nIndex,8,rs.strDeadRate);
}