大气扩散模型JAVA代码

package com; import java.util.ArrayList; import java.util.List; import com.mapinfo.mapj.MapJ; import com.mapinfo.util.DoublePoint; public class AirDensity { public double Q ; public double Pai = Math.PI; public double Gy ; public double Gz ; public double AirSpeed ; public double Dest; public double AirAngle ; public double X ; public double Y ; public double Z = 0;//; /** * @param q 单位排放量 * @param airSpeed 风速 * @param x 横坐标最大长度 * @param y 纵坐标最大长度 * @param airAngle 风向 */ public AirDensity(double q, double airSpeed, double x, double y,int airAngle) { Q = q; Dest=Math.sqrt(X * X + Y * Y); AirSpeed = airSpeed; AirAngle=airAngle*(Math.PI/180); //要转换成弧度 X = x; Y = y; } /**计算污染源浓度 * @return 浓度 */ public double CalDensity() { Gy = y1()*Math.pow(X,a1()) ; Gz = y2() * Math.pow(X, a2()); double density = Q * (1/(2 * Pai * AirSpeed*Gy*Gz))*Math.exp(-(Y*Y)/(2*(Gy*Gy)))*Math.exp(-0/(2*(Gz*Gz))); return density; } public double a1() { if (0 <= Dest && Dest <= 1000) return 0.926849; if ( Dest >1000) return 0.886940; return 0; } public double y1() { if (0 <= Dest && Dest <= 1000) return 0.143940; if (Dest > 1000) return 0.189396; return 0; } public double a2() { if (0 <= Dest && Dest <= 2000) return 0.838628; if (2000 < Dest && Dest <= 10000) return 0.756410; if (10000 < Dest ) return 0.815575; return 0; } public double y2() { if (0 <= Dest && Dest <= 2000) return 0.126152; if (2000 < Dest && Dest <= 10000) return 0.235667; if (10000 < Dest) return 0.136659; return 0; } /**根据x坐标和浓度计算y坐标 * @param x 横坐标 * @param density * @return */ public double CalYbyX(double x,double density) { AirDensity _Density=this;//new AirDensity(0,0,0,0,0); //_Density.Dest=x; double _a1=_Density.a1(); double _y1=_Density.y1(); double _a2=_Density.a2(); double _y2=_Density.y2(); double _Gy = _y1*Math.pow(x,_a1) ; double _Gz = _y2 * Math.pow(x, _a2); // double density = Q * (1/(2 * Pai * AirSpeed*Gy*Gz))*Math.exp(-(Y*Y)/(2*(Gy*Gy)))*Math.exp(-0/(2*(Gz*Gz))); double yy=density/Q; yy=yy*(2*Math.PI*AirSpeed*_Gy*_Gz); yy= Math.log(yy)/Math.log(Math.E) ; yy=Math.abs( yy*2*(Gy*Gy)); yy= Math.sqrt(yy); return yy; } private double miLon=0.0000115721704127876519013676624; //1米的经度 private double miLat=0.0000089932202939481227428007713; //1米的维度 /**按y轴步长200米获取10组点（每组11个) * @param orgPoint 原点坐标 * @return */ public List getPoints(DoublePoint orgPoint) { List<DensityPoints> densityPoints =new ArrayList<DensityPoints>(); for (int i = 1; i < 11; i++) { DensityPoints densityPoint =new DensityPoints(); DoublePoint[] dps=new DoublePoint[11]; int x=200*i; int y=0; this.X=x; this.Dest=x; double density=this.CalDensity(); densityPoint.setDensity(density); for(int j=0;j<6;j++){ double _x=x*(1-(1.0/15.0)*j); double _y=Math.abs( CalYbyX(_x,density)); double _cx1=_y; double _cy1=_x; //按风向旋转角度 //double cx1=_cx1*Math.cos(this.AirAngle) -_cy1*Math.sin(this.AirAngle); double cx1=_cy1*Math.sin(this.AirAngle)-_cx1*Math.cos(this.AirAngle) ; double cy1=_cx1*Math.sin(this.AirAngle)+_cy1*Math.cos(this.AirAngle); DoublePoint _orgPoint=new DoublePoint(orgPoint.x,orgPoint.y) ; _orgPoint.offset(cx1*miLon, cy1*miLat); dps[5-j]=new DoublePoint(_orgPoint.x,_orgPoint.y); if(j>0) { double _cx2=-_y; double _cy2=_x; // 按风向旋转角度 //double cx2=_cx2*Math.cos(this.AirAngle) -_cy2*Math.sin(this.AirAngle); double cx2=_cy2*Math.sin(this.AirAngle)-_cx2*Math.cos(this.AirAngle) ; double cy2=_cx2*Math.sin(this.AirAngle)+_cy2*Math.cos(this.AirAngle); _orgPoint=new DoublePoint(orgPoint.x,orgPoint.y) ; _orgPoint.offset(cx2*miLon, cy2*miLat); dps[5+j]=new DoublePoint(_orgPoint.x,_orgPoint.y); } } densityPoint.setPoints(dps); densityPoints.add(densityPoint); } return densityPoints; } }