fountain.cpp

此程序是用vc做的喷泉模拟实验 产生数据在砸爱matlab中画图就可以得到了

#include <iostream>
#include "Wind.h"
#include "Particle.h"
#include "Gravity.h"
#include <vector>
#include <time.h>
#include <fstream>
#include <math.h>
#include "cstdlib"
#define  num 1000
#define PI 3.14
//#define g 9.8

#include <string>
#include <sstream>

//#define PI 3.1415926  //pi值
//#define g 9.8 //重力加速度
#define num 200 //模拟水滴数目

using namespace std;

void main()
{   
	float V = 8.9; //喷头出水速度
	float vw = 1; //风速
	float aw = (float)90/180*PI; //风速与z轴夹角
	int count = 0; //落在花池外的水滴数据统计
	float R = 5;
	float r = 0.025;//喷水管半径
	float a = (float)15/180*PI; //喷出速度与Z轴的夹角
	CWind wind(vw,aw,PI/2);//pi/2定义为与x轴的夹角 
	wind.set_x(0);
	wind.set_y(vw*sin(aw));
	wind.set_z(vw*cos(aw));
	CGravity physics;
	std::vector<CParticle> drops;
	srand(time(0));
	for ( int i=0; i<num; i++ )
	{
		CParticle p(r*cos(2*PI/num*i), r*sin(2*PI/num*i),0.05);
		drops.push_back(p);
		float t = (float)rand()/RAND_MAX;
		drops.at(i).set_v((V+0.1*t)*sin(a)*cos(2*PI/num*i), (V+0.1*t)*sin(a)*sin(2*PI/num*i), (V+0.1*t)*cos(a));
		//设置水珠半径
		drops.at(i).set_r(0.0015+0.00005*t); 
	}
	float h0 = 0.05; //喷水管高度
	float ax = 0; float ay=0; float az=0; // 三个方向的加速度
	float step = 0.01;//模拟时间步长
	float c = 0.05;//水滴的空气阻力系数
	std::ofstream outfile("test.txt");
	std::vector<int> rem;
	rem.resize(100,0);
	while ( h0 > 0 )
	{  
		for ( int i=0; i<num; i++ )
		{
			//x轴方向上，只受阻力加速度影响，阻力加速度方向总是与物体运动速度方向相反
			if ( drops.at(i).get_vx() > 0 )
			{
				ax = -3*c*ro*pow((drops.at(i).get_vx()-wind.m_speedx),2)/4/drops.at(i).get_r();
			} 
			else if ( drops.at(i).get_vx() < 0 )
			{
				ax = 3*c*ro*pow((drops.at(i).get_vx()-wind.m_speedx),2)/4/drops.at(i).get_r();
			}
			else
			{
				ax = 0;
			}
			//y轴方向上，受阻力加速度和风力影响
			if ( drops.at(i).get_vy() > 0 )
			{
				ay = -3*c*ro*pow((drops.at(i).get_vy()-wind.m_speedy),2)/4/drops.at(i).get_r()-3*ro*pow(wind.m_speedy,2)/8/drops.at(i).get_r();
			}
			else if( drops.at(i).get_vy() < 0 )
			{
				ay = 3*c*ro*pow((drops.at(i).get_vy()-wind.m_speedy),2)/4/drops.at(i).get_r()-3*ro*pow(wind.m_speedy,2)/8/drops.at(i).get_r();
			}
			else
			{
				ay = -3*ro*pow(wind.m_speedy,2)/8/drops.at(i).get_r();
			}

		//z轴方向上，受风力和阻力的作用	

			if ( drops.at(i).get_vz() == 0)
			{
				az = -g-3*ro*pow(wind.m_speedz,2)/8/drops.at(i).get_r();
			} 
			else if ( drops.at(i).get_vz() > 0 )
			{
				az = -g-3*c*ro*pow(drops.at(i).get_vz()-wind.m_speedz,2)/8/drops.at(i).get_r()-3*ro*pow(wind.m_speedz,2)/8/drops.at(i).get_r();
			}
			else
			{
				az = -g+3*c*ro*pow(drops.at(i).get_vz()-wind.m_speedz,2)/8/drops.at(i).get_r()-3*ro*pow(wind.m_speedz,2)/8/drops.at(i).get_r();
			}


			float vx, vy, vz;
			float sx, sy, sz;
           
			vx = drops.at(i).get_vx()+ax*step;
			sx = drops.at(i).get_x()+drops.at(i).get_vx()*step+0.5*ax*pow(step,2);
			vy = drops.at(i).get_vy()+ay*step;
			sy = drops.at(i).get_y()+drops.at(i).get_vy()*step+0.5*ay*pow(step,2);

			vz = drops.at(i).get_vz()+az*step;
			sz = drops.at(i).get_z()+drops.at(i).get_vz()*step+0.5*az*pow(step,2);

			drops.at(i).set_position(sx, sy, sz);
			drops.at(i).set_v(vx, vy, vz);
			outfile<<sx<<" "<<sy<<"  "<<sz<<std::endl;
			if ( sqrt(pow(abs(sx),2)+pow(abs(sy),2)) > R && rem.at(i) == 0 )
			{
				rem.at(i) = 1; 
			}

			h0 = sz;
		}
		
	}

	for (  i=0; i<100; i++ )
	{
		if ( rem.at(i) == 1 )
		{
			count++;//统计落到花池外部的水滴个数 
		}
	}

	std::cout << count<<std::endl;


}