retrivelctrl.cpp

雷达资料的处理和应用

// RetrivelCtrl.cpp : implementation file
//

#include "stdafx.h"
#include "FuShe.h"
#include "RetrivelCtrl.h"
#include "Reference.h"
#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CRetrivelCtrl dialog


CRetrivelCtrl::CRetrivelCtrl(CWnd* pParent /*=NULL*/)
	: CDialog(CRetrivelCtrl::IDD, pParent)
{
	//{{AFX_DATA_INIT(CRetrivelCtrl)
		// NOTE: the ClassWizard will add member initialization here
	//}}AFX_DATA_INIT
}


void CRetrivelCtrl::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CRetrivelCtrl)
	DDX_Control(pDX, IDC_PROGRESS2, m_RetrieveCtrl);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CRetrivelCtrl, CDialog)
	//{{AFX_MSG_MAP(CRetrivelCtrl)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CRetrivelCtrl message handlers

BOOL CRetrivelCtrl::OnInitDialog() 
{
	CDialog::OnInitDialog();
	CProgressCtrl*Progress = (CProgressCtrl*)GetDlgItem(IDC_PROGRESS2);
	Progress->SetRange(1,nPoint*nPoint*4);
	Progress->SetStep(1);
	
	
	// TODO: Add extra initialization here
	
	return TRUE;  // return TRUE unless you set the focus to a control
	              // EXCEPTION: OCX Property Pages should return FALSE
}
BOOL CRetrivelCtrl::VelocityAzimuthProcess(CFile *file)
{
	unsigned char*ANGLE = new unsigned char[SWEEP*4];
	short int    *iBuff = new short int  [SWEEP*POINT];
	float        *fBuff = new float      [SWEEP*POINT];

	file->Read(m_HEAD,384);
	if ( m_HEAD[245] != 5 || m_HEAD[175] != 2)
		return FALSE;
	
	for (int k=0; k<SWEEP; k++)
	{
		file->Read(&ANGLE[k*4],4);
		m_fDIR[k] = (float)((ANGLE[k*4]*POINT
			+ANGLE[k*4+1])*Deg);
		m_fElevation += (float)((ANGLE[k*4+2]*POINT
			+ANGLE[k*4+3])*Deg);
		file->Read(&iBuff[k*POINT],POINT*2);
	}
	delete [] ANGLE;

	m_fElevation = (float)(m_fElevation/360.0-90.0);
	if (m_fElevation > 15.0) {
		AfxMessageBox("仰角过高");
		//return FALSE;
	}

	float*DEG = new float[SWEEP*POINT];
	float*VEL = new float[SWEEP*POINT];

	// The pulsation in the original data od Doppler velocity
	// has to smoothed before the retrieval.
	DimensionSmooth(iBuff,fBuff);
	delete [] iBuff;

	// Noted: calculate the horizonal velocity of 256 points on 
	// a sweep. the range_bin of doppler_radar is 0.25km and 
	// radius is 64.okm.
	for ( int j=0; j<POINT; j++ ) 
	for ( int i=0; i<SWEEP; i++ ) {// radar_sweep = 360
		int i1 = (SWEEP+i-1)%SWEEP;
		int i2 = (SWEEP+i+1)%SWEEP;
		if ( fBuff[i1*POINT+j]+fBuff[i2*POINT+j]==0.0 ) {
			DEG[i*POINT+j] = (float)0.0;
			VEL[i*POINT+j] = (float)0.0;
		}
		else {
			DEG[i*POINT+j] = (float)(atan(
				-(fBuff[i1*POINT+j]-fBuff[i2*POINT+j])
				/((fBuff[i1*POINT+j]+fBuff[i2*POINT+j])
				*tan(PI))));
			if (fabs(DEG[i*POINT+j]/PI)<85.0) {
				VEL[i*POINT+j] = (float)( LSB*fabs(
				      (fBuff[i1*POINT+j]+fBuff[i2*POINT+j])
					  /(2*cos(DEG[i*POINT+j])*cos(PI))));
			}
			else
				VEL[i*POINT+j] = (float)0.0;
		}
		if ( fBuff[i1*POINT+j]+fBuff[i2*POINT+j]<0 ) {
			m_fUwind[i*POINT+j] = (float)
				( VEL[i*POINT+j]*sin(DEG[i*POINT+j]+m_fDIR[i]*PI) );
			m_fVwind[i*POINT+j] = (float)
				( VEL[i*POINT+j]*cos(DEG[i*POINT+j]+m_fDIR[i]*PI) );
		}
		else {
			m_fUwind[i*POINT+j] = (float)
				( VEL[i*POINT+j]*
				  sin(DEG[i*POINT+j]+m_fDIR[i]*PI+PAI) );
			m_fVwind[i*POINT+j] = (float)
				( VEL[i*POINT+j]*
				  cos(DEG[i*POINT+j]+m_fDIR[i]*PI+PAI) );
		}
	}
	delete [] DEG;
	delete [] VEL;

	delete [] fBuff;

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// The pulsation in the original data od Doppler velocity has to 
// smoothed before the retrieval. But it can be brought into being by
// one-demension linear smooth.
BOOL CRetrivelCtrl::DimensionSmooth(short int*iBuff,float*fBuff)
{
	for (int Point=0; Point<POINT; Point++) {
		for (int i=0; i<SWEEP; i++) {
			int i1 = (SWEEP+i-1)%SWEEP;
			int i2 = (SWEEP+i+1)%SWEEP;

			fBuff[i*POINT+Point] = (float)(
				(1.0-Weigh)*iBuff[i*POINT+Point]
				+(iBuff[i2*POINT+Point]
				+iBuff[i1*POINT+Point])
				*(Weigh/2.0));
		}
	}

	return TRUE;
}

/////////////////////////////////////////////////////////////////////////
// cx,cy:   分别是直角坐标下的横坐标和纵坐标;
// radius:  是cressman半径,作为插值的标准;
// distance:两点间的距离;
// weight:  权重系数;
BOOL CRetrivelCtrl::CressmanInterval()
{
	float radius,weight,distance,cx[nPoint*2],cy[nPoint*2];

	int *AX = new int[SWEEP*POINT];
	int *AY = new int[SWEEP*POINT];
	GetPoint(AX,AY);
	
	for(int i=0; i<nPoint; ++i) {
		cx[i] = (float)(scale/nPoint*(i-nPoint));
		cy[i] = (float)(scale/nPoint*(nPoint-i));
        
		cx[nPoint+i] = (float)(scale/nPoint*i);
        cy[nPoint+i] =-(float)(scale/nPoint*i);
    }
    radius = (float)(scale/nPoint*Factor);
    for( i=0; i<4*nPoint*nPoint; ++i ) {
		m_RetrieveCtrl.SetPos(i);
		int x=i%(2*nPoint);
		int y=(i-i%(2*nPoint))/(2*nPoint);
		m_Uwind[i] = (float)0.0;
        m_Vwind[i] = (float)0.0;
        weight  = (float)0.0;
        for(int k=0; k<SWEEP*POINT; ++k) {
            distance = (float)( sqrt((AX[k]-cx[x])*(AX[k]-cx[x])
				+(AY[k]-cy[y])*(AY[k]-cy[y])));
            if( distance < radius ) {
                m_Uwind[i]+=(float)((radius-distance)/(radius+distance)*m_fUwind[k]);
				m_Vwind[i]+=(float)((radius-distance)/(radius+distance)*m_fVwind[k]);
                weight+=(radius-distance)/(radius+distance);
            }
        }
		m_Uwind[i]=(float)(m_Uwind[i]/weight);
        m_Vwind[i]=(float)(m_Vwind[i]/weight);
	}

	m_RetrieveCtrl.SetPos(0);

	delete [] AX;
	delete [] AY;

	return TRUE;
}

BOOL CRetrivelCtrl::GetPoint(int *X, int *Y)
{
	double  SIN,   COS;

	for (int i=0; i<SWEEP; i++) {
		SIN=sin(m_fDIR[i]*PI);
		COS=cos(m_fDIR[i]*PI);
		for (int j=0; j<POINT; j++) {
			X[i*POINT+j]=(int)(SIN*j*0.25);
			Y[i*POINT+j]=(int)(COS*j*0.25);
		}
	}

	return TRUE;
}

BOOL CRetrivelCtrl::DrawWindField(CClientDC*pdc)
{
	int x_begin = 100, y_begin = 100;

	if (DrawBox(x_begin,y_begin,pdc) == FALSE)
		return FALSE;

	for (int Point=0; Point<nPoint*nPoint*4; Point++) {
		int x = Point%(2*nPoint);
		int y = (Point-Point%(2*nPoint))/(2*nPoint);
		
		if ((x>0 && x<2*nPoint-1) && (y>0 && y<2*nPoint-1)) {
			int   x_tail,y_tail,x_head,y_head;
			float direction,uu,vv,windLength;
			y_tail = y_begin+y*yLength;
			x_tail = x_begin+x*xLength;
			pdc->SetPixel(x_tail,y_tail,RGB(255,255,255));
		
			uu = m_Uwind[Point];
			vv = m_Vwind[Point];//????
			windLength = (float)(sqrt(uu*uu+vv*vv));
			if (windLength >1.0&&windLength <40.0) {
				direction = ReturnDirection(uu,vv);
				x_head = x_tail+(int)(sin(direction)*windLength
					*ArrowLength/UnitLength);
				y_head = y_tail-(int)(cos(direction)*windLength
					*ArrowLength/UnitLength);
			
				// draw the length of the wind vector:
				pdc->MoveTo(x_tail,y_tail);
				pdc->LineTo(x_head,y_head);

				// draw the tail of the wind vector:
				DrawWindTail(pdc,x_head,y_head,direction);
			}
		}
	}

	return TRUE;
}

float CRetrivelCtrl::ReturnDirection(float u, float v)
{
	float  direction = (float)0.0;

	if (v != 0.0) {
		if ( v > 0.0 ) 
			direction = (float)(2*PAI+atan(u/v));
		else
			direction = (float)(PAI+atan(u/v));
	}
	else {
		if ( u > 0.0 )
			direction = (float)(PAI/2.0);
		if ( u < 0.0 )
			direction = (float)(3*PAI/2.0);
	}

	return direction;
}

BOOL CRetrivelCtrl::DrawBox(int x, int y, CClientDC*pdc)
{
	if (x<0 || x>800 || y<0 || y>600)
		return FALSE;

	int x_end = x+(2*nPoint)*xLength,y_end = y+(2*nPoint)*yLength;
	pdc->MoveTo(x,y);
	pdc->LineTo(x_end,y);
	pdc->TextOut(x-50,y,"60.0");

	pdc->MoveTo(x_end,y);
	pdc->LineTo(x_end,y_end);
	pdc->TextOut(x_end,y_end+15,"54.0");

	pdc->MoveTo(x_end,y_end);
	pdc->LineTo(x,y_end);
	pdc->TextOut(x-15,y_end+15,"-60.0");

	pdc->MoveTo(x,y_end);
	pdc->LineTo(x,y);
	pdc->TextOut(x-50,y_end-5,"-54.0");

	for (int i=0; i<2*nPoint; i++) {
		int x1, y1, x2, y2;
		x1 = x;
		y1 = (int)(y+i*yLength);
		x2 = x-8;
		y2 = (int)(y+i*yLength);
		pdc->MoveTo(x1,y1);
		pdc->LineTo(x2,y2);

		y1 = y_end;
		x1 = (int)(x+i*xLength);
		y2 = y_end+8;
		x2 = (int)(x+i*xLength);
		pdc->MoveTo(x1,y1);
		pdc->LineTo(x2,y2);
	}
	pdc->TextOut(x-30,y+yLength*nPoint," 0");
	pdc->TextOut(x+xLength*nPoint-5,y_end+15," 0");
	pdc->TextOut((x+x_end)/2,y_end+50,"(km)");
	// draw unit:
	pdc->MoveTo(x+125,y-10);
	pdc->LineTo(x+125,y-15);
	pdc->MoveTo(x+125,y-10);
	pdc->LineTo(x+125+ArrowLength,y-10);
	pdc->MoveTo(x+125+ArrowLength,y-10);
	pdc->LineTo(x+125+ArrowLength,y-15);
	pdc->TextOut(x+150,y-20,"15(m/s)");
	
	// display time, date and elevation
	char datetime[20],csElevation[10];
	sprintf(datetime,"%02d\\%02d\\%02d    %02d:%02d:%02d",m_HEAD[63],
		m_HEAD[67],m_HEAD[71],m_HEAD[75],m_HEAD[79],m_HEAD[83]);
	pdc->TextOut(x+30,y-60,datetime);
	sprintf(csElevation,"%5.2f",m_fElevation);
	pdc->TextOut(x+200,y-60,"Ele:");
	pdc->TextOut(x+220,y-60,csElevation);

	return TRUE;
}

BOOL CRetrivelCtrl::DrawWindTail(CClientDC*pdc,int x,int y,float dir)
{
	double dArrowAngle = 0.42;
	int x_arrow, y_arrow;

	x_arrow = x+(int)(-sin(dir+dArrowAngle)*5.0);
    y_arrow = y+(int)(cos(dir+dArrowAngle)*5.0);
	
	pdc->MoveTo(x,y);
	pdc->LineTo(x_arrow,y_arrow);

	x_arrow = x+(int)(-sin(dir-dArrowAngle)*5.0);
    y_arrow = y+(int)(cos(dir-dArrowAngle)*5.0);

	pdc->MoveTo(x,y);
	pdc->LineTo(x_arrow,y_arrow);

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// There are two method(VAd and VAP) to retrieve the wind vetor field 
// based on single-doppler velocity. 
// VAD(Velocity Azimuth Display) method is based on the assumption the
// horizonal wind vector are homogeneity. 
BOOL CRetrivelCtrl::VelocityAzimuthDisplay(CFile*file,int*cn)
{
	unsigned char*ANGLE = new unsigned char  [SWEEP*4];
	short int    *iBuff = new short int  [SWEEP*POINT];
	float        *fBuff = new float      [SWEEP*POINT];

	file->Read(m_HEAD,384);
	if ( m_HEAD[245] != 5 || m_HEAD[175] != 2)
		return FALSE;
	
	for (int k=0; k<SWEEP; k++)
	{
		file->Read(&ANGLE[k*4],4);
		m_fDIR[k] = (float)((ANGLE[k*4]*POINT
			+ANGLE[k*4+1])*Deg);
		m_fElevation += (float)((ANGLE[k*4+2]*POINT
			+ANGLE[k*4+3])*Deg);
		file->Read(&iBuff[k*POINT],POINT*2);
	}
	delete [] ANGLE;

	m_fElevation = (float)(m_fElevation/360.0-90.0);

	// The pulsation in the original data od Doppler velocity
	// has to smoothed before the retrieval.
	DimensionSmooth(iBuff,fBuff);
	delete [] iBuff;

	// Noted: calculate the horizonal velocity of 256 points on 
	// a sweep. the range_bin of doppler_radar is 0.25km and 
	// radius is 64.okm.
	for ( int j=0; j<POINT; j++ ) 
	{
		double tempU = 0.0, tempV = 0.0;
		for ( int i=0; i<SWEEP/2; i++ ) // radar_sweep = 360
		{
			tempU += fBuff[i*POINT+j]*sin(m_fDIR[i]*PI)*PI;
			tempV += fBuff[i*POINT+j]*cos(m_fDIR[i]*PI)*PI;
		}
		m_fUwind[*cn*POINT+j] = (float)(tempU*LSB);
		m_fVwind[*cn*POINT+j] = (float)(tempV*LSB);
	}
	delete [] fBuff;

	return TRUE;
}