喷泉模拟仿真的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 100 //模拟水滴数目 using namespace std; void main() { float V = 8.9; //喷头出水速度 float vw = 0; //风速 float aw = (float)90/180*PI; //风速与z轴夹角 int count = 0; //落在花池外的水滴数据统计 float R = 2.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; }