一些常用的物理量计算

气象领域中一些常用的物理量计算代码，如水汽压、饱和水汽压、混合比、位温、假相当位温、比湿、饱和比湿、露点温度、水的汽化潜热、相对湿度、绝对湿度等。以及球面距离公式。

---

Formula.java

---

001  /**
002     PACKAGE  : cma.common.atmos
003 
004          功能：气象领域的常用物理量计算
005        文件名：Formula.java
006      相关文件：
007        编写者：刘泽军
008      开始日期：
009      更新日期：
010          备注：
011 
012    */
013 
014  package  cma.common.atmos;
015 
016  import  java.io.*;
017  import  java.lang.*;
018  import  java.util.*;
019 
020  public class  Formula  {
021  /**
022    * 构造函数
023    */
024       public  Formula () {
025       }
026       public static  double  VALUE  =  9999 ; //Algorithm.DefaultValue;    //缺省值
027  /**
028    * 功能：计算两点距离
029    * 参数：
030    * lon1, lat1   - 第1点位置(度)
031    * lon2, lat2   - 第2点位置(度)
032    * 返回值：
033    *      两点距离
034    */
035       public static  double  distanceDegree ( double  lon1,  double  lat1,  double  lon2,  double  lat2 ) {
036           return ( distance ( Math.toRadians ( lon1 ) , Math.toRadians ( lat1 ) , Math.toRadians ( lon1 ) , Math.toRadians ( lat1 ))) ;
037       }
038  /**
039    * 功能：计算两点距离
040    * 参数：
041    * lon1, lat1   - 第1点位置(弧度)
042    * lon2, lat2   - 第2点位置(弧度)
043    * 返回值：
044    *      两点距离
045    */
046       public static  double  distance ( double  lon1,  double  lat1,  double  lon2,  double  lat2 ) { //(弧度)
047           return ( Meteorology.RADIUS * Math.acos ( Math.sin ( lat1 ) *Math.sin ( lat2 )  + Math.cos ( lat1 ) *Math.cos ( lat2 ) *Math.cos ( lon2-lon1 ))) ;
048           //地球半径*arccsos(sin(wd1)sin(wd2)+cos(wd1)cos(wd2)cos(jd2-jd1)) 球面距离？
049       }
050       public static  double  w ( double  _p /*气压*/ ,  double  _e /*水汽压*/ ) { //混合比
051           if (  VALUE != _p && VALUE != _e &&  0  != _p && _p != _e  ) {
052               return ( 0.62198  * _e /  (  _p - _e  )) ;
053           }
054           return ( VALUE ) ;
055       }
056       public static  double  e ( double  _p /*气压*/ ,  double  _q /*比湿*/ ) { //水汽压
057           if (  VALUE != _p && VALUE != _q &&  0  != _p  ) {
058               return ( _p * _q /  0.62198 ) ;
059           }
060           return ( VALUE ) ;
061       }
062       public static  double  e ( double  _td /*露点温度*/ ) { //水汽压
063           //《大气科学常用公式》P382 td => e
064           if (  _td ==  0.0  ) {
065               return ( 6.1078 ) ;
066           }
067           else if (  VALUE != _td && - 237.3  != _td  ) {
068               return ( 6.1078  * Math.pow ( 10 ,  7.5 / ( 237.3 /_td+ 1 ))) ;
069           }
070           else  {
071               return ( VALUE ) ;
072           }
073       }
074       public static  double  E ( double  _t /*温度*/ ) { //饱和水汽压
075           //《大气科学常用公式》P2
076           if (  VALUE != _t && - 243.92  != _t  ) {
077               return ( 6.1078  * Math.pow ( 10 , 7.69 *_t/ ( _t+ 243.92 ))) ;
078           }
079           return ( VALUE ) ;
080       }
081       public static  double  a ( double  _e /*水汽压*/ ,  double  _t /*温度*/ ) { //绝对湿度=e/Rv/T
082           double   _Rv  =  461.51 ; //J . kg-1 . K-1
083           if (  VALUE != _e && VALUE != _t && - 273.16  != _t  ) {
084               return ( _e /  ( _Rv* ( 273.16 +_t ))) ;
085           }
086           return ( VALUE ) ;
087       }
088       public static  double  r ( double  _e /*水汽压*/ ,  double  _E /*饱和水汽压*/ ) { //相对湿度=e/E*100%
089           if (  VALUE != _e && VALUE != _E  ) {
090               return ( _e / _E *  100.0 ) ;
091           }
092           return ( VALUE ) ;
093       }
094       public static  double  theta ( double  _p /*气压*/ ,  double  _t /*温度*/ ) { //理想气体的位温
095           //以干空气的位温代替
096           return ( thetad ( _p, _t )) ;
097       }
098       public static  double  thetad ( double  _p /*气压*/ ,  double  _t /*温度*/ ) { //干空气的位温
099           if (  VALUE != _t && VALUE != _p &&  0.0  != _p  ) {
100               return ( ( _t +  273.16 )  * Math.pow ( 1000 /_p,  0.286 ) ) ; //《大气科学常用公式》P466 pd= p
101           }
102           return ( VALUE ) ;
103       }
104       public static  double  thetam ( double  _p /*气压*/ ,  double  _t /*温度*/ ,  double  _e /*水汽压*/ ,  double  _q /*比湿*/ ) { //湿空气的位温
105           double   _Rm      =  287.04  *  (  1  +  0.608  * _q  ) ;
106           double   _Cpm     =  1005  *  (  1  +  0.86  * _q  ) ;
107           if (  VALUE != _t && VALUE != _p &&  0.0  != _p && VALUE != _e &&  0.0  != _p + _e && VALUE != _q  ) {
108               return ( ( _t +  273.16 )  * Math.pow ( 1000 / ( _p+_e ) ,_Rm/_Cpm ) ) ; //《大气科学常用公式》P466
109           }
110           return ( VALUE ) ;
111       }
112       public static  double  blowupTheta1 ( double  _p /*气压*/ ,  double  _t /*温度*/ ,  double  _td /*露点温度*/ ) { //溃变 - θ1
113           return ( theta ( _p, _t )) ;
114       }
115       public static  double  blowupTheta2 ( double  _p /*气压*/ ,  double  _t /*温度*/ ,  double  _td /*露点温度*/ ) { //溃变 - θ2
116           if (  _p ==  0  || _p == VALUE || _t == VALUE || _td == VALUE  )   return ( VALUE ) ;
117           double   _o1 = blowupTheta1 ( _p, _t, _td ) ;
118           if (  VALUE == _o1  ) {
119               return ( VALUE ) ;
120           }
121           double   _h  =  124  *  (  _t - _td  ) ;
122           double   _tc = _t -  0.0098  * _h;
123           double   _e  =  6.1078  * Math.pow (  10 ,  ( 7.5 *_tc ) / ( 237.3 +_tc ) ) ;
124           double   _qs =  0.62198  * _e / _p;
125           double   _cp =  1005 ;
126           double   _l  =  2500000 ;
127           double   _o2 = _o1 * Math.exp ( ( _l*_qs )  /  ( _cp* ( _tc+ 273.16 )) ) ;
128           return ( _o2 ) ;
129       }
130       public static  double   blowupTheta3 ( double  _p /*气压*/ ,  double  _t /*温度*/ ,  double  _td /*露点温度*/ ) { //溃变 - θ3
131           if (  _p ==  0  || _p == VALUE || _t == VALUE || _td == VALUE  )   return ( VALUE ) ;
132           double   _o1 = blowupTheta1 ( _p, _t, _td ) ;
133           if (  _o1 == VALUE  )  return ( VALUE ) ;
134           double   _e1 =  6.1078  * Math.pow (  10 ,  ( 7.5 *_t ) / ( 237.3 +_t ) ) ;
135           double   _cp =  1005 ;
136           double   _q  =  0.62198  * _e1 / _p;
137           double   _l  =  2500000 ;
138           double   _o3 = _o1 * Math.exp (( _l*_q ) / ( _cp* ( _t+ 273.16 ))) ;
139           return (  _o3  ) ;
140       }
141       public static  double   thetase ( double  _p /*气压*/ ,  double  _t /*温度*/ ) { //假相当位温
142           double   thetase = VALUE;
143           return ( thetase ) ;
144       }
145       public static  double   q ( double  _e /*水汽压*/ ,  double  _p /*气压*/ ) { //比湿(e为饱和水汽压返回饱和比湿，否则返回比湿)
146           if (  VALUE != _e && VALUE != _p &&  0.0  != _p  ) {
147               return ( 0.62198  * _e / _p ) ; //《大气科学常用公式》P382
148           }
149           return ( VALUE ) ;
150       }
151       public static  double  qs ( double  _E /*饱和水汽压*/ ,  double  _p /*气压*/ ) { //饱和比湿
152           return ( q ( _E, _p )) ;
153       }
154       public static  double  dv ( double  _dd /*风向*/ ,  double  _ff /*风速*/ ) { //风的垂直分量
155           double   vv = VALUE;                                              //    ^
156           if (  VALUE != _dd && VALUE != _ff  ) {                             //   /|/N
157                                                                           //W   |   E
158               vv = Math.abs (  _ff * Math.cos ( Math.toRadians ( _dd ))) ;         //----+--->
159               if (  _dd >=  270.0  || _dd <  90.0  ) {                           //    |
160                   vv = - 1.0  * Math.abs ( vv ) ;                                //    | S
161               }                                                            //    |
162           }
163           return ( vv ) ;
164       }
165       public static  double  du ( double  _dd /*风向*/ ,  double  _ff /*风速*/ ) { //风的水平分量
166           double   _vh = VALUE;                                             //    ^
167           if (  VALUE != _dd && VALUE != _ff  ) {                             //   /|/N
168                                                                           //W   |   E
169               _vh = Math.abs (  _ff * Math.sin ( Math.toRadians ( _dd )) ) ;       //----+--->
170               if (  _dd <=  180.0  ) {                                         //    |
171                   _vh = - Math.abs (  _vh  ) ;                                 //    | S
172               }                                                            //    |
173           }
174           return ( _vh ) ;
175       }
176       public static  double   td ( double  _e /*水汽压*/ ) { //利用水汽压计算露点温度
177           double   _td = VALUE;
178           double   _A  =  6.1078 , _a =  7.5 , _b =  237.3 ;
179           if (  _e == VALUE  ) {
180               _td = VALUE;
181           }
182           else if (  _e == _A  ) {
183               _td =  0 ;
184           }
185           else  {
186               _td = _b /  (  _a / Math.log ( _e/_A ) /Math.log ( 10.0 )  -  1 ) ; //《大气科学常用公式》P382
187           }
188           return ( _td ) ;
189       }
190       public static  double   td ( double  _p /*气压*/ ,double _t /*温度*/ ,double _r /*相对湿度*/ ) { //利用气压、气温及相对湿度计算露点温度 - 待检查
191           double   _E   = E ( _t ) ; // t => E
192           double   _e   = _E == VALUE ? VALUE : _E * _r /  100.0 ; // 《大气科学常用公式》P381 r = e / E * 100 %
193           return (  td ( _e ) ) ;
194       }
195       public static  double   lv ( double  _t /*温度*/ ) //水的汽化潜热
196  { //单位：J . kg -1
197           double   _Rv  =  461.51 ; //J . kg-1 . K -1
198           double   _a   =  7.726 ,    _a1  =  7.44 ;
199           double   _b   =  244.88 ,   _b1  =  235.33 ;
200           double   _lv  = VALUE;
201           if (  _t != VALUE && _t <=  0  && _t != -_b  ) {
202               _lv = _a * _b * _Rv * Math.log ( 10 )  * Math.pow ( ( _t +  237.16 ) / ( _b+_t ) ,  2  ) ;
203           }
204           else if (  _t != VALUE && _t >  0  && _t != -_b1  ) {
205               _lv = _a1 * _b1 * _Rv * Math.log ( 10 )  * Math.pow ( ( _t +  237.16 ) / ( _b1+_t ) ,  2  ) ;
206           }
207           return ( _lv ) ;
208       }
209  /*
210       public static String LongLat2dfm(double value) {
211           double  d1,f1,m1;
212           TStringList *StringList1 = new TStringList();
213           TStringList *StringList2 = new TStringList();
214           FountainExtractString(FormatFloat("0.000000",value),".",false,StringList1);
215           FountainExtractString(FormatFloat("0.000000",(AnsiString("0.")+StringList1->Strings[1]).ToDouble()*60),".",false,StringList2);
216 
217           d1 = StringList1->Strings[0].ToDouble();
218           f1 = StringList2->Strings[0].ToDouble();
219           m1 = (AnsiString("0.")+StringList2->Strings[1]).ToDouble()*60;
220 
221           if( m1 >= 59.5 ) {
222               m1 = 0.0;
223               f1 += 1.0;
224           }
225           if( f1 >= 59.5 ) {
226               f1 = 0.0;
227               d1 += 1.0;
228           }
229 
230           AnsiString FReturnStr =
231                   FormatFloat("0",d1)+"°"+
232                   FormatFloat("00",f1)+"′"+
233                   FormatFloat("00",m1)+"″";
234           delete StringList1;
235           delete StringList2;
236           return(FReturnStr);
237       }
238  */
239  }