サンプルJAVAプログラム( メタン (CH4) )

メタンのページに戻る
トップページ(2電子原子も含む新ボーア模型)

このサンプルプログラムでは、実行すると、まず、電子a (ea) の最初のx座標 r (MM) を入力するように画面に表示される。
電子b(水素原子の)に関しては、最初の座標 aa と bb (in MM) を入力する。
aa は 水素原子核 b と point B (=電子 b の座標を 炭素原子核と水素原子核 b を結ぶ線に投影した交点) の間の距離。
(aa が負の値のとき、この point B と 原点は 水素原子核 b を挟んで互いに逆に位置する。)
bb は、この point B と 電子 b の間の距離。(この bb が正の値のとき、電子 b のy座標は、負になる。)

この入力値から、炭素原子の4つの外殻軌道の電子 (ea, ec, ef,eh) と、4つの水素原子の電子 (eb, ed, ee, eg) が正四面体状に配置される。
ここでは、便宜のため、新しい単位を使っている。( 1 MM = 10-14 meter, 1 SS = 10-21 second, 1 MM/SS = 107 m/s )
電子a が 自身の1.0, 1.9, 2.0 ド・ブロイ波長分の軌道を進み、かつ、電子 b が 1ド・ブロイ波長分の軌道を進んだとき、このプログラムは次の値を画面に表示する。
各電子の位置座標、
最初と最後の地点における炭素(もしくは水素)原子核と 電子a(もしくは電子b) の距離、
また、電子a の軌道上の点の中で、最も核との距離が長く(もしくは短く)なるときの、その距離。


import java.util.Scanner;
class MathMethod {
 public static void main(String[] args) {
 
 Scanner stdIn=new Scanner(System.in);     
 System.out.println("initial x coordinate r of electron a (MM)? ");  
 double r=stdIn.nextDouble();
 System.out.println("aa= eb distance 1 from H+ nucleus (MM)? ");  
 double aa=stdIn.nextDouble();
 System.out.println("bb = eb distance 2 from H+ nucleus (MM)? ");  
 double bb=stdIn.nextDouble();
 
                         // E=CH4 ionization energy
 double E=219.6753; double EE=-E*1.60217646e-19;
 double me=9.1093826e-31; double n=4;    // n=number of carbon valence electrons
 double pai=3.141592653589793; double epsi=8.85418781787346e-12;
 double h=6.62606896e-34; double ele=1.60217653e-19;
 double rsi=Math.sqrt(6);double rtw=Math.sqrt(2);double rth=Math.sqrt(3);
 
 double ka,kb,kc,kd,ke,kf,kg,kh;
 double ahb,ahd,ahe,ahg,bhb,bhd,bhe,bhg,chb,chd,che,chg,dhb,dhd,dhe,dhg;
 double ehb,ehd,ehe,ehg,fhb,fhd,fhe,fhg,ghb,ghd,ghe,ghg,hhb,hhd,hhe,hhg;
 double hbx,hby,hbz,hdx,hdy,hdz,hex,hey,hez,hgx,hgy,hgz;
 double kab,kac,kad,kae,kaf,kag,kah,kbc,kbd,kbe,kbf,kbg,kbh,kcd,kce,kcf,kcg,kch,kde,kdf,kdg,kdh,kef,keg,keh,kfg,kfh,kgh;
 double ax,ay,az,pax,pay,paz,vax,vay,vaz,wa,vka;
 double bx,by,bz,pbx,pby,pbz,vbx,vby,vbz,wb,vkb;
 double cx,cy,cz,pcx,pcy,pcz,vcx,vcy,vcz,wc,vkc;
 double dx,dy,dz,pdx,pdy,pdz,vdx,vdy,vdz,wd,vkd;
 double ex,ey,ez,pex,pey,pez,vex,vey,vez,we,vke;
 double fx,fy,fz,pfx,pfy,pfz,vfx,vfy,vfz,wf,vkf;
 double gx,gy,gz,pgx,pgy,pgz,vgx,vgy,vgz,wg,vkg;
 double hx,hy,hz,phx,phy,phz,vhx,vhy,vhz,wh,vkh;
    
 double rr=r*1.0e-14; double bbb=(rsi*r)/2;     // bbb=initial distance between nucleus and electron
 double sz=r/rtw;                               // sz=initial z-coordinate of electron a
                                         // hydrogen atom (hb,hd,he,hg) positions 
 hbx=8899.8;hby=0.0;hbz=6293.1;hdx=-hbx;hdy=-hby;hdz=hbz;hex=hby;hey=-hbx;hez=-hbz;
 hgx=-hby;hgy=hbx;hgz=-hbz;
 
 double kra=Math.sqrt(hbx*hbx+hby*hby+hbz*hbz);double krb=Math.sqrt(hbx*hbx+hby*hby);
                          // set initial positions of electrons (b,d,e,g) from inputted aa and bb
 pbx= hbx-(aa*hbx)/kra+(hby*bb)/krb; pby=hby-(hby*aa)/kra-(hbx*bb)/krb; pbz=hbz-(hbz*aa)/kra;
 pdx= hdx-(aa*hdx)/kra+(hdy*bb)/krb; pdy=hdy-(hdy*aa)/kra-(hdx*bb)/krb; pdz=hdz-(hdz*aa)/kra;
 pex= hex-(aa*hex)/kra+(hey*bb)/krb; pey=hey-(hey*aa)/kra-(hex*bb)/krb; pez=hez-(hez*aa)/kra;
 pgx= hgx-(aa*hgx)/kra+(hgy*bb)/krb; pgy=hgy-(hgy*aa)/kra-(hgx*bb)/krb; pgz=hgz-(hgz*aa)/kra;
                                    // show initial coordinate of electron b
 System.out.printf("initial bx:%.2f ", pbx);System.out.printf("initial by:%.2f ", pby);System.out.printf("initial bz:%.2f \n", pbz);
                                   // set initial positions of electrons (a, c, f, h)
 pax=r;pay=0.0;paz=sz; pcx=-r;pcy=0.0;pcz=sz; 
 pfx=0.0;pfy=-r;pfz=-sz; phx=0.0;phy=r;phz=-sz;
                                  // show initial coordinate of electron a
 System.out.printf("initial ax:%.2f ", pax);System.out.printf("initial ay:%.2f ", pay);System.out.printf("initial az:%.2f \n", paz);
                                 // show initial distance between electron b and hydrogen b
 double hyko=Math.sqrt((hbx-pbx)*(hbx-pbx)+(hby-pby)*(hby-pby)+(hbz-pbz)*(hbz-pbz));
 System.out.printf("initial distance between eb and Hb:%.2f \n", hyko);
 
                                // distances between electrons and hydrogen nuclei
 ahb=Math.sqrt((pax-hbx)*(pax-hbx)+(pay-hby)*(pay-hby)+(paz-hbz)*(paz-hbz));
 ahd=Math.sqrt((pax-hdx)*(pax-hdx)+(pay-hdy)*(pay-hdy)+(paz-hdz)*(paz-hdz));
 ahe=Math.sqrt((pax-hex)*(pax-hex)+(pay-hey)*(pay-hey)+(paz-hez)*(paz-hez));
 ahg=Math.sqrt((pax-hgx)*(pax-hgx)+(pay-hgy)*(pay-hgy)+(paz-hgz)*(paz-hgz));
 bhb=Math.sqrt((pbx-hbx)*(pbx-hbx)+(pby-hby)*(pby-hby)+(pbz-hbz)*(pbz-hbz));
 bhd=Math.sqrt((pbx-hdx)*(pbx-hdx)+(pby-hdy)*(pby-hdy)+(pbz-hdz)*(pbz-hdz));
 bhe=Math.sqrt((pbx-hex)*(pbx-hex)+(pby-hey)*(pby-hey)+(pbz-hez)*(pbz-hez));
 bhg=Math.sqrt((pbx-hgx)*(pbx-hgx)+(pby-hgy)*(pby-hgy)+(pbz-hgz)*(pbz-hgz));
 chb=Math.sqrt((pcx-hbx)*(pcx-hbx)+(pcy-hby)*(pcy-hby)+(pcz-hbz)*(pcz-hbz));
 chd=Math.sqrt((pcx-hdx)*(pcx-hdx)+(pcy-hdy)*(pcy-hdy)+(pcz-hdz)*(pcz-hdz));
 che=Math.sqrt((pcx-hex)*(pcx-hex)+(pcy-hey)*(pcy-hey)+(pcz-hez)*(pcz-hez));
 chg=Math.sqrt((pcx-hgx)*(pcx-hgx)+(pcy-hgy)*(pcy-hgy)+(pcz-hgz)*(pcz-hgz));
 dhb=Math.sqrt((pdx-hbx)*(pdx-hbx)+(pdy-hby)*(pdy-hby)+(pdz-hbz)*(pdz-hbz));
 dhd=Math.sqrt((pdx-hdx)*(pdx-hdx)+(pdy-hdy)*(pdy-hdy)+(pdz-hdz)*(pdz-hdz));
 dhe=Math.sqrt((pdx-hex)*(pdx-hex)+(pdy-hey)*(pdy-hey)+(pdz-hez)*(pdz-hez));
 dhg=Math.sqrt((pdx-hgx)*(pdx-hgx)+(pdy-hgy)*(pdy-hgy)+(pdz-hgz)*(pdz-hgz));
 ehb=Math.sqrt((pex-hbx)*(pex-hbx)+(pey-hby)*(pey-hby)+(pez-hbz)*(pez-hbz));
 ehd=Math.sqrt((pex-hdx)*(pex-hdx)+(pey-hdy)*(pey-hdy)+(pez-hdz)*(pez-hdz));
 ehe=Math.sqrt((pex-hex)*(pex-hex)+(pey-hey)*(pey-hey)+(pez-hez)*(pez-hez));
 ehg=Math.sqrt((pex-hgx)*(pex-hgx)+(pey-hgy)*(pey-hgy)+(pez-hgz)*(pez-hgz));
 fhb=Math.sqrt((pfx-hbx)*(pfx-hbx)+(pfy-hby)*(pfy-hby)+(pfz-hbz)*(pfz-hbz));
 fhd=Math.sqrt((pfx-hdx)*(pfx-hdx)+(pfy-hdy)*(pfy-hdy)+(pfz-hdz)*(pfz-hdz));
 fhe=Math.sqrt((pfx-hex)*(pfx-hex)+(pfy-hey)*(pfy-hey)+(pfz-hez)*(pfz-hez));
 fhg=Math.sqrt((pfx-hgx)*(pfx-hgx)+(pfy-hgy)*(pfy-hgy)+(pfz-hgz)*(pfz-hgz));
 ghb=Math.sqrt((pgx-hbx)*(pgx-hbx)+(pgy-hby)*(pgy-hby)+(pgz-hbz)*(pgz-hbz));
 ghd=Math.sqrt((pgx-hdx)*(pgx-hdx)+(pgy-hdy)*(pgy-hdy)+(pgz-hdz)*(pgz-hdz));
 ghe=Math.sqrt((pgx-hex)*(pgx-hex)+(pgy-hey)*(pgy-hey)+(pgz-hez)*(pgz-hez));
 ghg=Math.sqrt((pgx-hgx)*(pgx-hgx)+(pgy-hgy)*(pgy-hgy)+(pgz-hgz)*(pgz-hgz));
 hhb=Math.sqrt((phx-hbx)*(phx-hbx)+(phy-hby)*(phy-hby)+(phz-hbz)*(phz-hbz));
 hhd=Math.sqrt((phx-hdx)*(phx-hdx)+(phy-hdy)*(phy-hdy)+(phz-hdz)*(phz-hdz));
 hhe=Math.sqrt((phx-hex)*(phx-hex)+(phy-hey)*(phy-hey)+(phz-hez)*(phz-hez));
 hhg=Math.sqrt((phx-hgx)*(phx-hgx)+(phy-hgy)*(phy-hgy)+(phz-hgz)*(phz-hgz));

                          // distance between each electron and carbon nucleus
 ka=Math.sqrt(pax*pax+pay*pay+paz*paz);
 kb=Math.sqrt(pbx*pbx+pby*pby+pbz*pbz);
 kc=Math.sqrt(pcx*pcx+pcy*pcy+pcz*pcz);
 kd=Math.sqrt(pdx*pdx+pdy*pdy+pdz*pdz);
 ke=Math.sqrt(pex*pex+pey*pey+pez*pez);
 kf=Math.sqrt(pfx*pfx+pfy*pfy+pfz*pfz);
 kg=Math.sqrt(pgx*pgx+pgy*pgy+pgz*pgz);
 kh=Math.sqrt(phx*phx+phy*phy+phz*phz);

                           // distances between two electrons
 kab=Math.sqrt((pax-pbx)*(pax-pbx)+(pay-pby)*(pay-pby)+(paz-pbz)*(paz-pbz));
 kac=Math.sqrt((pax-pcx)*(pax-pcx)+(pay-pcy)*(pay-pcy)+(paz-pcz)*(paz-pcz));
 kad=Math.sqrt((pax-pdx)*(pax-pdx)+(pay-pdy)*(pay-pdy)+(paz-pdz)*(paz-pdz));
 kae=Math.sqrt((pax-pex)*(pax-pex)+(pay-pey)*(pay-pey)+(paz-pez)*(paz-pez)); 
 kaf=Math.sqrt((pax-pfx)*(pax-pfx)+(pay-pfy)*(pay-pfy)+(paz-pfz)*(paz-pfz));
 kag=Math.sqrt((pax-pgx)*(pax-pgx)+(pay-pgy)*(pay-pgy)+(paz-pgz)*(paz-pgz));
 kah=Math.sqrt((pax-phx)*(pax-phx)+(pay-phy)*(pay-phy)+(paz-phz)*(paz-phz));
 kbc=Math.sqrt((pbx-pcx)*(pbx-pcx)+(pby-pcy)*(pby-pcy)+(pbz-pcz)*(pbz-pcz)); 
 kbd=Math.sqrt((pbx-pdx)*(pbx-pdx)+(pby-pdy)*(pby-pdy)+(pbz-pdz)*(pbz-pdz));
 kbe=Math.sqrt((pbx-pex)*(pbx-pex)+(pby-pey)*(pby-pey)+(pbz-pez)*(pbz-pez));
 kbf=Math.sqrt((pbx-pfx)*(pbx-pfx)+(pby-pfy)*(pby-pfy)+(pbz-pfz)*(pbz-pfz));
 kbg=Math.sqrt((pbx-pgx)*(pbx-pgx)+(pby-pgy)*(pby-pgy)+(pbz-pgz)*(pbz-pgz));
 kbh=Math.sqrt((pbx-phx)*(pbx-phx)+(pby-phy)*(pby-phy)+(pbz-phz)*(pbz-phz));
 kcd=Math.sqrt((pcx-pdx)*(pcx-pdx)+(pcy-pdy)*(pcy-pdy)+(pcz-pdz)*(pcz-pdz));
 kce=Math.sqrt((pcx-pex)*(pcx-pex)+(pcy-pey)*(pcy-pey)+(pcz-pez)*(pcz-pez));
 kcf=Math.sqrt((pcx-pfx)*(pcx-pfx)+(pcy-pfy)*(pcy-pfy)+(pcz-pfz)*(pcz-pfz));
 kcg=Math.sqrt((pcx-pgx)*(pcx-pgx)+(pcy-pgy)*(pcy-pgy)+(pcz-pgz)*(pcz-pgz));
 kch=Math.sqrt((pcx-phx)*(pcx-phx)+(pcy-phy)*(pcy-phy)+(pcz-phz)*(pcz-phz));
 kde=Math.sqrt((pdx-pex)*(pdx-pex)+(pdy-pey)*(pdy-pey)+(pdz-pez)*(pdz-pez));
 kdf=Math.sqrt((pdx-pfx)*(pdx-pfx)+(pdy-pfy)*(pdy-pfy)+(pdz-pfz)*(pdz-pfz));
 kdg=Math.sqrt((pdx-pgx)*(pdx-pgx)+(pdy-pgy)*(pdy-pgy)+(pdz-pgz)*(pdz-pgz));
 kdh=Math.sqrt((pdx-phx)*(pdx-phx)+(pdy-phy)*(pdy-phy)+(pdz-phz)*(pdz-phz));
 kef=Math.sqrt((pex-pfx)*(pex-pfx)+(pey-pfy)*(pey-pfy)+(pez-pfz)*(pez-pfz));  
 keg=Math.sqrt((pex-pgx)*(pex-pgx)+(pey-pgy)*(pey-pgy)+(pez-pgz)*(pez-pgz));
 keh=Math.sqrt((pex-phx)*(pex-phx)+(pey-phy)*(pey-phy)+(pez-phz)*(pez-phz));
 kfg=Math.sqrt((pfx-pgx)*(pfx-pgx)+(pfy-pgy)*(pfy-pgy)+(pfz-pgz)*(pfz-pgz)); 
 kfh=Math.sqrt((pfx-phx)*(pfx-phx)+(pfy-phy)*(pfy-phy)+(pfz-phz)*(pfz-phz));  
 kgh=Math.sqrt((pgx-phx)*(pgx-phx)+(pgy-phy)*(pgy-phy)+(pgz-phz)*(pgz-phz));

                          // pta,ptb= potential energy of electron a,b
 double ac=(ele*ele)/(4*pai*epsi);
 double pta=ac*(-n/ka+1/kab+1/kac+1/kad+1/kae+1/kaf+1/kag+1/kah-1/ahb-1/ahd-1/ahe-1/ahg);
 double ptb=ac*(-n/kb+1/kab+1/kbc+1/kbd+1/kbe+1/kbf+1/kbg+1/kbh-1/bhb-1/bhd-1/bhe-1/bhg);
 
 double krd=Math.sqrt((hbx-hdx)*(hbx-hdx)+(hby-hdy)*(hby-hdy)+(hbz-hdz)*(hbz-hdz));
                           // poten=total potential energy
 double poten=ac*(-n/ka-n/kb-n/kc-n/kd-n/ke-n/kf-n/kg-n/kh+1/kab+1/kac+1/kad+1/kae+1/kaf+1/kag+1/kah+
1/kbc+1/kbd+1/kbe+1/kbf+1/kbg+1/kbh+1/kcd+1/kce+1/kcf+1/kcg+1/kch+1/kde
+1/kdf+1/kdg+1/kdh+1/kef+1/keg+1/keh+1/kfg+1/kfh+1/kgh-1/ahb-1/ahd-1/ahe-1/ahg-1/bhb-1/bhd-1/bhe-1/bhg-1/chb-1/chd
-1/che-1/chg-1/dhb-1/dhd-1/dhe-1/dhg-1/ehb-1/ehd-1/ehe-1/ehg-1/fhb-1/fhd-1/fhe-1/fhg-1/ghb-1/ghd-1/ghe-1/ghg
-1/hhb-1/hhd-1/hhe-1/hhg+(4*n)/kra+6/krd); 

 poten=poten*1.0e14; pta=pta*1.0e14;ptb=ptb*1.0e14;
                         
 double TT=EE-poten;        // distribute  kinetic energy based on each potential energy 
 double Ta=(TT*pta)/(4*(pta+ptb)); double Tb=(TT*ptb)/(4*(pta+ptb)); 

                                 // change MM to meter
 hbx=hbx*1.0e-14;hby=hby*1.0e-14;hbz=hbz*1.0e-14;
 hdx=hdx*1.0e-14;hdy=hdy*1.0e-14;hdz=hdz*1.0e-14;
 hex=hex*1.0e-14;hey=hey*1.0e-14;hez=hez*1.0e-14;
 hgx=hgx*1.0e-14;hgy=hgy*1.0e-14;hgz=hgz*1.0e-14;

 if (Ta>0 && Tb>0) {
 double vya=Math.sqrt((2*Ta)/(me));
 double vv=vya*1.0e-7;         // change m/s to MM/SS
                               // set each initial velocity vector from kinetic energy
 double hita=Math.sqrt(65);   
 vax=0.0;vaz=-vv/hita;vay=(vv*8)/hita; vcx=-vax; vcy=-vay; vcz=vaz; vfx=vay; vfy=-vax; vfz=-vaz;
  vhx=-vay; vhy=vax; vhz=-vaz;
 double vyb=Math.sqrt((2*Tb)/(me));
 double vvb=vyb*1.0e-7; 
 
 vbz=vvb; vbx=0.0; vby=0.0;
 vez=-vvb; vex=0.0; vey=0.0;
 vgz=-vvb; vgx=0.0; vgy=0.0;
 vdz=vvb; vdx=0.0; vdy=0.0; 
       
 wa=0.0;wb=0.0;wc=0.0;wd=0.0;we=0.0;wf=0.0;wg=0.0;wh=0.0; ax=pax; ay=pay;az=paz;
 double bigdis=Math.sqrt(ax*ax+ay*ay+az*az);double shortdis=bigdis; double kko; 
  double hyou=-1.0;double hhaa=-1.0;  double hrr=-1.0;   // mark
  double hrrr=-1.0;double hhrr=-1.0;double saigoh;
  
 do {
    ax=pax+vax; ay=pay+vay; az=paz+vaz; bx=pbx+vbx; by=pby+vby; bz=pbz+vbz;cx=pcx+vcx; cy=pcy+vcy; cz=pcz+vcz;
    dx=pdx+vdx; dy=pdy+vdy; dz=pdz+vdz; ex=pex+vex; ey=pey+vey; ez=pez+vez; fx=pfx+vfx; fy=pfy+vfy; fz=pfz+vfz;
    gx=pgx+vgx; gy=pgy+vgy; gz=pgz+vgz; hx=phx+vhx; hy=phy+vhy; hz=phz+vhz;
                                   
    kko=Math.sqrt(pax*pax+pay*pay+paz*paz);
 
    if (bigdis < kko) {bigdis=kko;}

    if (shortdis > kko) {shortdis=kko;} 

                                      // show each electron's position on the way of the orbit
    if (wb>1.0 && hrr<0) {hrr=1.0; System.out.printf("wb:%.2f \n", wb);System.out.printf("ax:%.2f ", ax);System.out.printf("ay:%.2f ", ay);System.out.printf("az:%.2f \n", az);
  System.out.printf("bx:%.2f ", bx);System.out.printf("by:%.2f ", by);System.out.printf("bz:%.2f \n", bz);
  System.out.printf("cx:%.2f ", cx);System.out.printf("cy:%.2f ", cy);System.out.printf("cz:%.2f \n", cz);
  System.out.printf("dx:%.2f ", dx);System.out.printf("dy:%.2f ", dy);System.out.printf("dz:%.2f \n", dz);
  System.out.printf("ex:%.2f ", ex);System.out.printf("ey:%.2f ", ey);System.out.printf("ez:%.2f \n", ez);
  System.out.printf("fx:%.2f ", fx);System.out.printf("fy:%.2f ", fy);System.out.printf("fz:%.2f \n", fz);
  System.out.printf("gx:%.2f ", gx);System.out.printf("gy:%.2f ", gy);System.out.printf("gz:%.2f \n", gz);
  System.out.printf("hx:%.2f ", hx);System.out.printf("hy:%.2f ", hy);System.out.printf("hz:%.2f \n", hz);
   saigoh=Math.sqrt((hbx*1.0e14-pbx)*(hbx*1.0e14-pbx)+(hby*1.0e14-pby)*(hby*1.0e14-pby)+(hbz*1.0e14-pbz)*(hbz*1.0e14-pbz));
 System.out.printf("last distance between Hb and eb:%.2f\n ", saigoh); }

  if (wa>2.0 && hrrr<0) {hrrr=1.0; 
   System.out.printf("wa:%.2f \n", wa);
  System.out.printf("ax:%.2f ", ax);System.out.printf("ay:%.2f ", ay);System.out.printf("az:%.2f \n", az);
  System.out.printf("bx:%.2f ", bx);System.out.printf("by:%.2f ", by);System.out.printf("bz:%.2f \n", bz);
  System.out.printf("cx:%.2f ", cx);System.out.printf("cy:%.2f ", cy);System.out.printf("cz:%.2f \n", cz);
  System.out.printf("dx:%.2f ", dx);System.out.printf("dy:%.2f ", dy);System.out.printf("dz:%.2f \n", dz);
  System.out.printf("ex:%.2f ", ex);System.out.printf("ey:%.2f ", ey);System.out.printf("ez:%.2f \n", ez);
  System.out.printf("fx:%.2f ", fx);System.out.printf("fy:%.2f ", fy);System.out.printf("fz:%.2f \n", fz);
  System.out.printf("gx:%.2f ", gx);System.out.printf("gy:%.2f ", gy);System.out.printf("gz:%.2f \n", gz);
  System.out.printf("hx:%.2f ", hx);System.out.printf("hy:%.2f ", hy);System.out.printf("hz:%.2f \n", hz);
  System.out.printf("initial distance between nucleus and ea:%.2f \n ", bbb);
  System.out.printf("last distance between nucleus and ea :%.2f  \n  ",   kko);}

 if (hrr>0 && hrrr>0) {hhrr=1.0;}

    if (wa>1.0 && hyou<0) {hyou=1.0; System.out.printf("wa:%.2f \n", wa);System.out.printf("ax:%.2f ", ax);System.out.printf("ay:%.2f ", ay);System.out.printf("az:%.2f \n", az);
  System.out.printf("bx:%.2f ", bx);System.out.printf("by:%.2f ", by);System.out.printf("bz:%.2f \n", bz);
  System.out.printf("cx:%.2f ", cx);System.out.printf("cy:%.2f ", cy);System.out.printf("cz:%.2f \n", cz);
  System.out.printf("dx:%.2f ", dx);System.out.printf("dy:%.2f ", dy);System.out.printf("dz:%.2f \n", dz);
  System.out.printf("ex:%.2f ", ex);System.out.printf("ey:%.2f ", ey);System.out.printf("ez:%.2f \n", ez);
  System.out.printf("fx:%.2f ", fx);System.out.printf("fy:%.2f ", fy);System.out.printf("fz:%.2f \n", fz);
  System.out.printf("gx:%.2f ", gx);System.out.printf("gy:%.2f ", gy);System.out.printf("gz:%.2f \n", gz);
  System.out.printf("hx:%.2f ", hx);System.out.printf("hy:%.2f ", hy);System.out.printf("hz:%.2f \n", hz);}


    if (wa>1.9 && hhaa<0) {hhaa=1.0; System.out.printf("wa:%.2f \n", wa);System.out.printf("ax:%.2f ", ax);System.out.printf("ay:%.2f ", ay);System.out.printf("az:%.2f \n", az);
  System.out.printf("bx:%.2f ", bx);System.out.printf("by:%.2f ", by);System.out.printf("bz:%.2f \n", bz);
  System.out.printf("cx:%.2f ", cx);System.out.printf("cy:%.2f ", cy);System.out.printf("cz:%.2f \n", cz);
  System.out.printf("dx:%.2f ", dx);System.out.printf("dy:%.2f ", dy);System.out.printf("dz:%.2f \n", dz);
  System.out.printf("ex:%.2f ", ex);System.out.printf("ey:%.2f ", ey);System.out.printf("ez:%.2f \n", ez);
  System.out.printf("fx:%.2f ", fx);System.out.printf("fy:%.2f ", fy);System.out.printf("fz:%.2f \n", fz);
  System.out.printf("gx:%.2f ", gx);System.out.printf("gy:%.2f ", gy);System.out.printf("gz:%.2f \n", gz);
  System.out.printf("hx:%.2f ", hx);System.out.printf("hy:%.2f ", hy);System.out.printf("hz:%.2f \n", hz);}

                               // calculate de Broglie's waves
    vka=vax*vax+vay*vay+vaz*vaz; wa=wa+(me*vka*1.0e-7)/h; vkb=vbx*vbx+vby*vby+vbz*vbz; wb=wb+(me*vkb*1.0e-7)/h;
    vkc=vcx*vcx+vcy*vcy+vcz*vcz; wc=wc+(me*vkc*1.0e-7)/h; vkd=vdx*vdx+vdy*vdy+vdz*vdz; wd=wd+(me*vkd*1.0e-7)/h;
    vke=vex*vex+vey*vey+vez*vez; we=we+(me*vke*1.0e-7)/h; vkf=vfx*vfx+vfy*vfy+vfz*vfz; wf=wf+(me*vkf*1.0e-7)/h; 
    vkg=vgx*vgx+vgy*vgy+vgz*vgz; wg=wg+(me*vkg*1.0e-7)/h; vkh=vhx*vhx+vhy*vhy+vhz*vhz; wh=wh+(me*vkh*1.0e-7)/h;
                                  
                               // change MM to meter
    pax=pax*1.0e-14;pay=pay*1.0e-14;paz=paz*1.0e-14; pbx=pbx*1.0e-14;pby=pby*1.0e-14;pbz=pbz*1.0e-14;
    pcx=pcx*1.0e-14;pcy=pcy*1.0e-14;pcz=pcz*1.0e-14; pdx=pdx*1.0e-14;pdy=pdy*1.0e-14;pdz=pdz*1.0e-14;
    pex=pex*1.0e-14;pey=pey*1.0e-14;pez=pez*1.0e-14; pfx=pfx*1.0e-14;pfy=pfy*1.0e-14;pfz=pfz*1.0e-14;
    pgx=pgx*1.0e-14;pgy=pgy*1.0e-14;pgz=pgz*1.0e-14; phx=phx*1.0e-14;phy=phy*1.0e-14;phz=phz*1.0e-14;

                               // distance between carbon nucleus and each electron
    ka=Math.sqrt(pax*pax+pay*pay+paz*paz);kb=Math.sqrt(pbx*pbx+pby*pby+pbz*pbz);
    kc=Math.sqrt(pcx*pcx+pcy*pcy+pcz*pcz);kd=Math.sqrt(pdx*pdx+pdy*pdy+pdz*pdz);
    ke=Math.sqrt(pex*pex+pey*pey+pez*pez);kf=Math.sqrt(pfx*pfx+pfy*pfy+pfz*pfz);
    kg=Math.sqrt(pgx*pgx+pgy*pgy+pgz*pgz);kh=Math.sqrt(phx*phx+phy*phy+phz*phz);

                               // distances between electrons and hydrogen nuclei
 ahb=Math.sqrt((pax-hbx)*(pax-hbx)+(pay-hby)*(pay-hby)+(paz-hbz)*(paz-hbz));
 ahd=Math.sqrt((pax-hdx)*(pax-hdx)+(pay-hdy)*(pay-hdy)+(paz-hdz)*(paz-hdz));
 ahe=Math.sqrt((pax-hex)*(pax-hex)+(pay-hey)*(pay-hey)+(paz-hez)*(paz-hez));
 ahg=Math.sqrt((pax-hgx)*(pax-hgx)+(pay-hgy)*(pay-hgy)+(paz-hgz)*(paz-hgz));
 bhb=Math.sqrt((pbx-hbx)*(pbx-hbx)+(pby-hby)*(pby-hby)+(pbz-hbz)*(pbz-hbz));
 bhd=Math.sqrt((pbx-hdx)*(pbx-hdx)+(pby-hdy)*(pby-hdy)+(pbz-hdz)*(pbz-hdz));
 bhe=Math.sqrt((pbx-hex)*(pbx-hex)+(pby-hey)*(pby-hey)+(pbz-hez)*(pbz-hez));
 bhg=Math.sqrt((pbx-hgx)*(pbx-hgx)+(pby-hgy)*(pby-hgy)+(pbz-hgz)*(pbz-hgz));
 chb=Math.sqrt((pcx-hbx)*(pcx-hbx)+(pcy-hby)*(pcy-hby)+(pcz-hbz)*(pcz-hbz));
 chd=Math.sqrt((pcx-hdx)*(pcx-hdx)+(pcy-hdy)*(pcy-hdy)+(pcz-hdz)*(pcz-hdz));
 che=Math.sqrt((pcx-hex)*(pcx-hex)+(pcy-hey)*(pcy-hey)+(pcz-hez)*(pcz-hez));
 chg=Math.sqrt((pcx-hgx)*(pcx-hgx)+(pcy-hgy)*(pcy-hgy)+(pcz-hgz)*(pcz-hgz));
 dhb=Math.sqrt((pdx-hbx)*(pdx-hbx)+(pdy-hby)*(pdy-hby)+(pdz-hbz)*(pdz-hbz));
 dhd=Math.sqrt((pdx-hdx)*(pdx-hdx)+(pdy-hdy)*(pdy-hdy)+(pdz-hdz)*(pdz-hdz));
 dhe=Math.sqrt((pdx-hex)*(pdx-hex)+(pdy-hey)*(pdy-hey)+(pdz-hez)*(pdz-hez));
 dhg=Math.sqrt((pdx-hgx)*(pdx-hgx)+(pdy-hgy)*(pdy-hgy)+(pdz-hgz)*(pdz-hgz));
 ehb=Math.sqrt((pex-hbx)*(pex-hbx)+(pey-hby)*(pey-hby)+(pez-hbz)*(pez-hbz));
 ehd=Math.sqrt((pex-hdx)*(pex-hdx)+(pey-hdy)*(pey-hdy)+(pez-hdz)*(pez-hdz));
 ehe=Math.sqrt((pex-hex)*(pex-hex)+(pey-hey)*(pey-hey)+(pez-hez)*(pez-hez));
 ehg=Math.sqrt((pex-hgx)*(pex-hgx)+(pey-hgy)*(pey-hgy)+(pez-hgz)*(pez-hgz));
 fhb=Math.sqrt((pfx-hbx)*(pfx-hbx)+(pfy-hby)*(pfy-hby)+(pfz-hbz)*(pfz-hbz));
 fhd=Math.sqrt((pfx-hdx)*(pfx-hdx)+(pfy-hdy)*(pfy-hdy)+(pfz-hdz)*(pfz-hdz));
 fhe=Math.sqrt((pfx-hex)*(pfx-hex)+(pfy-hey)*(pfy-hey)+(pfz-hez)*(pfz-hez));
 fhg=Math.sqrt((pfx-hgx)*(pfx-hgx)+(pfy-hgy)*(pfy-hgy)+(pfz-hgz)*(pfz-hgz));
 ghb=Math.sqrt((pgx-hbx)*(pgx-hbx)+(pgy-hby)*(pgy-hby)+(pgz-hbz)*(pgz-hbz));
 ghd=Math.sqrt((pgx-hdx)*(pgx-hdx)+(pgy-hdy)*(pgy-hdy)+(pgz-hdz)*(pgz-hdz));
 ghe=Math.sqrt((pgx-hex)*(pgx-hex)+(pgy-hey)*(pgy-hey)+(pgz-hez)*(pgz-hez));
 ghg=Math.sqrt((pgx-hgx)*(pgx-hgx)+(pgy-hgy)*(pgy-hgy)+(pgz-hgz)*(pgz-hgz));
 hhb=Math.sqrt((phx-hbx)*(phx-hbx)+(phy-hby)*(phy-hby)+(phz-hbz)*(phz-hbz));
 hhd=Math.sqrt((phx-hdx)*(phx-hdx)+(phy-hdy)*(phy-hdy)+(phz-hdz)*(phz-hdz));
 hhe=Math.sqrt((phx-hex)*(phx-hex)+(phy-hey)*(phy-hey)+(phz-hez)*(phz-hez));
 hhg=Math.sqrt((phx-hgx)*(phx-hgx)+(phy-hgy)*(phy-hgy)+(phz-hgz)*(phz-hgz));

                                      // distances between two electrons
 kab=Math.sqrt((pax-pbx)*(pax-pbx)+(pay-pby)*(pay-pby)+(paz-pbz)*(paz-pbz));
 kac=Math.sqrt((pax-pcx)*(pax-pcx)+(pay-pcy)*(pay-pcy)+(paz-pcz)*(paz-pcz));
 kad=Math.sqrt((pax-pdx)*(pax-pdx)+(pay-pdy)*(pay-pdy)+(paz-pdz)*(paz-pdz));
 kae=Math.sqrt((pax-pex)*(pax-pex)+(pay-pey)*(pay-pey)+(paz-pez)*(paz-pez)); 
 kaf=Math.sqrt((pax-pfx)*(pax-pfx)+(pay-pfy)*(pay-pfy)+(paz-pfz)*(paz-pfz));
 kag=Math.sqrt((pax-pgx)*(pax-pgx)+(pay-pgy)*(pay-pgy)+(paz-pgz)*(paz-pgz));
 kah=Math.sqrt((pax-phx)*(pax-phx)+(pay-phy)*(pay-phy)+(paz-phz)*(paz-phz));
 kbc=Math.sqrt((pbx-pcx)*(pbx-pcx)+(pby-pcy)*(pby-pcy)+(pbz-pcz)*(pbz-pcz)); 
 kbd=Math.sqrt((pbx-pdx)*(pbx-pdx)+(pby-pdy)*(pby-pdy)+(pbz-pdz)*(pbz-pdz));
 kbe=Math.sqrt((pbx-pex)*(pbx-pex)+(pby-pey)*(pby-pey)+(pbz-pez)*(pbz-pez));
 kbf=Math.sqrt((pbx-pfx)*(pbx-pfx)+(pby-pfy)*(pby-pfy)+(pbz-pfz)*(pbz-pfz));
 kbg=Math.sqrt((pbx-pgx)*(pbx-pgx)+(pby-pgy)*(pby-pgy)+(pbz-pgz)*(pbz-pgz));
 kbh=Math.sqrt((pbx-phx)*(pbx-phx)+(pby-phy)*(pby-phy)+(pbz-phz)*(pbz-phz));
 kcd=Math.sqrt((pcx-pdx)*(pcx-pdx)+(pcy-pdy)*(pcy-pdy)+(pcz-pdz)*(pcz-pdz));
 kce=Math.sqrt((pcx-pex)*(pcx-pex)+(pcy-pey)*(pcy-pey)+(pcz-pez)*(pcz-pez));
 kcf=Math.sqrt((pcx-pfx)*(pcx-pfx)+(pcy-pfy)*(pcy-pfy)+(pcz-pfz)*(pcz-pfz));
 kcg=Math.sqrt((pcx-pgx)*(pcx-pgx)+(pcy-pgy)*(pcy-pgy)+(pcz-pgz)*(pcz-pgz));
 kch=Math.sqrt((pcx-phx)*(pcx-phx)+(pcy-phy)*(pcy-phy)+(pcz-phz)*(pcz-phz));
 kde=Math.sqrt((pdx-pex)*(pdx-pex)+(pdy-pey)*(pdy-pey)+(pdz-pez)*(pdz-pez));
 kdf=Math.sqrt((pdx-pfx)*(pdx-pfx)+(pdy-pfy)*(pdy-pfy)+(pdz-pfz)*(pdz-pfz));
 kdg=Math.sqrt((pdx-pgx)*(pdx-pgx)+(pdy-pgy)*(pdy-pgy)+(pdz-pgz)*(pdz-pgz));
 kdh=Math.sqrt((pdx-phx)*(pdx-phx)+(pdy-phy)*(pdy-phy)+(pdz-phz)*(pdz-phz));
 kef=Math.sqrt((pex-pfx)*(pex-pfx)+(pey-pfy)*(pey-pfy)+(pez-pfz)*(pez-pfz));  
 keg=Math.sqrt((pex-pgx)*(pex-pgx)+(pey-pgy)*(pey-pgy)+(pez-pgz)*(pez-pgz));
 keh=Math.sqrt((pex-phx)*(pex-phx)+(pey-phy)*(pey-phy)+(pez-phz)*(pez-phz));
 kfg=Math.sqrt((pfx-pgx)*(pfx-pgx)+(pfy-pgy)*(pfy-pgy)+(pfz-pgz)*(pfz-pgz)); 
 kfh=Math.sqrt((pfx-phx)*(pfx-phx)+(pfy-phy)*(pfy-phy)+(pfz-phz)*(pfz-phz));  
 kgh=Math.sqrt((pgx-phx)*(pgx-phx)+(pgy-phy)*(pgy-phy)+(pgz-phz)*(pgz-phz));
                             
                                     //  change each velocity vector 
     double ai=(ele*ele)/(4*pai*epsi*me);
    vax=vax+1.0e-28*ai*(-(n*pax)/(ka*ka*ka)+(pax-pbx)/(kab*kab*kab)+(pax-pcx)/(kac*kac*kac)+(pax-pdx)/(kad*kad*kad)+(pax-pex)/(kae*kae*kae)+(pax-pfx)/(kaf*kaf*kaf)+(pax-pgx)/(kag*kag*kag)+(pax-phx)/(kah*kah*kah)-(pax-hbx)/(ahb*ahb*ahb)-(pax-hdx)/(ahd*ahd*ahd)-(pax-hex)/(ahe*ahe*ahe)-(pax-hgx)/(ahg*ahg*ahg));  
     vay=vay+1.0e-28*ai*(-(n*pay)/(ka*ka*ka)+(pay-pby)/(kab*kab*kab)+(pay-pcy)/(kac*kac*kac)+(pay-pdy)/(kad*kad*kad)+(pay-pey)/(kae*kae*kae)+(pay-pfy)/(kaf*kaf*kaf)+(pay-pgy)/(kag*kag*kag)+(pay-phy)/(kah*kah*kah)-(pay-hby)/(ahb*ahb*ahb)-(pay-hdy)/(ahd*ahd*ahd)-(pay-hey)/(ahe*ahe*ahe)-(pay-hgy)/(ahg*ahg*ahg));
     vaz=vaz+1.0e-28*ai*(-(n*paz)/(ka*ka*ka)+(paz-pbz)/(kab*kab*kab)+(paz-pcz)/(kac*kac*kac)+(paz-pdz)/(kad*kad*kad)+(paz-pez)/(kae*kae*kae)+(paz-pfz)/(kaf*kaf*kaf)+(paz-pgz)/(kag*kag*kag)+(paz-phz)/(kah*kah*kah)-(paz-hbz)/(ahb*ahb*ahb)-(paz-hdz)/(ahd*ahd*ahd)-(paz-hez)/(ahe*ahe*ahe)-(paz-hgz)/(ahg*ahg*ahg)); 
     
    vbx=vbx+1.0e-28*ai*(-(n*pbx)/(kb*kb*kb)+(pbx-pax)/(kab*kab*kab)+(pbx-pcx)/(kbc*kbc*kbc)+(pbx-pdx)/(kbd*kbd*kbd)+(pbx-pex)/(kbe*kbe*kbe)+(pbx-pfx)/(kbf*kbf*kbf)+(pbx-pgx)/(kbg*kbg*kbg)+(pbx-phx)/(kbh*kbh*kbh)-(pbx-hbx)/(bhb*bhb*bhb)-(pbx-hdx)/(bhd*bhd*bhd)-(pbx-hex)/(bhe*bhe*bhe)-(pbx-hgx)/(bhg*bhg*bhg));  
     vby=vby+1.0e-28*ai*(-(n*pby)/(kb*kb*kb)+(pby-pay)/(kab*kab*kab)+(pby-pcy)/(kbc*kbc*kbc)+(pby-pdy)/(kbd*kbd*kbd)+(pby-pey)/(kbe*kbe*kbe)+(pby-pfy)/(kbf*kbf*kbf)+(pby-pgy)/(kbg*kbg*kbg)+(pby-phy)/(kbh*kbh*kbh)-(pby-hby)/(bhb*bhb*bhb)-(pby-hdy)/(bhd*bhd*bhd)-(pby-hey)/(bhe*bhe*bhe)-(pby-hgy)/(bhg*bhg*bhg));
     vbz=vbz+1.0e-28*ai*(-(n*pbz)/(kb*kb*kb)+(pbz-paz)/(kab*kab*kab)+(pbz-pcz)/(kbc*kbc*kbc)+(pbz-pdz)/(kbd*kbd*kbd)+(pbz-pez)/(kbe*kbe*kbe)+(pbz-pfz)/(kbf*kbf*kbf)+(pbz-pgz)/(kbg*kbg*kbg)+(pbz-phz)/(kbh*kbh*kbh)-(pbz-hbz)/(bhb*bhb*bhb)-(pbz-hdz)/(bhd*bhd*bhd)-(pbz-hez)/(bhe*bhe*bhe)-(pbz-hgz)/(bhg*bhg*bhg)); 

   vcx=vcx+1.0e-28*ai*(-(n*pcx)/(kc*kc*kc)+(pcx-pax)/(kac*kac*kac)+(pcx-pbx)/(kbc*kbc*kbc)+(pcx-pdx)/(kcd*kcd*kcd)+(pcx-pex)/(kce*kce*kce)+(pcx-pfx)/(kcf*kcf*kcf)+(pcx-pgx)/(kcg*kcg*kcg)+(pcx-phx)/(kch*kch*kch)-(pcx-hbx)/(chb*chb*chb)-(pcx-hdx)/(chd*chd*chd)-(pcx-hex)/(che*che*che)-(pcx-hgx)/(chg*chg*chg));  
     vcy=vcy+1.0e-28*ai*(-(n*pcy)/(kc*kc*kc)+(pcy-pay)/(kac*kac*kac)+(pcy-pby)/(kbc*kbc*kbc)+(pcy-pdy)/(kcd*kcd*kcd)+(pcy-pey)/(kce*kce*kce)+(pcy-pfy)/(kcf*kcf*kcf)+(pcy-pgy)/(kcg*kcg*kcg)+(pcy-phy)/(kch*kch*kch)-(pcy-hby)/(chb*chb*chb)-(pcy-hdy)/(chd*chd*chd)-(pcy-hey)/(che*che*che)-(pcy-hgy)/(chg*chg*chg));
     vcz=vcz+1.0e-28*ai*(-(n*pcz)/(kc*kc*kc)+(pcz-paz)/(kac*kac*kac)+(pcz-pbz)/(kbc*kbc*kbc)+(pcz-pdz)/(kcd*kcd*kcd)+(pcz-pez)/(kce*kce*kce)+(pcz-pfz)/(kcf*kcf*kcf)+(pcz-pgz)/(kcg*kcg*kcg)+(pcz-phz)/(kch*kch*kch)-(pcz-hbz)/(chb*chb*chb)-(pcz-hdz)/(chd*chd*chd)-(pcz-hez)/(che*che*che)-(pcz-hgz)/(chg*chg*chg));   

   vdx=vdx+1.0e-28*ai*(-(n*pdx)/(kd*kd*kd)+(pdx-pax)/(kad*kad*kad)+(pdx-pbx)/(kbd*kbd*kbd)+(pdx-pcx)/(kcd*kcd*kcd)+(pdx-pex)/(kde*kde*kde)+(pdx-pfx)/(kdf*kdf*kdf)+(pdx-pgx)/(kdg*kdg*kdg)+(pdx-phx)/(kdh*kdh*kdh)-(pdx-hbx)/(dhb*dhb*dhb)-(pdx-hdx)/(dhd*dhd*dhd)-(pdx-hex)/(dhe*dhe*dhe)-(pdx-hgx)/(dhg*dhg*dhg));  
     vdy=vdy+1.0e-28*ai*(-(n*pdy)/(kd*kd*kd)+(pdy-pay)/(kad*kad*kad)+(pdy-pby)/(kbd*kbd*kbd)+(pdy-pcy)/(kcd*kcd*kcd)+(pdy-pey)/(kde*kde*kde)+(pdy-pfy)/(kdf*kdf*kdf)+(pdy-pgy)/(kdg*kdg*kdg)+(pdy-phy)/(kdh*kdh*kdh)-(pdy-hby)/(dhb*dhb*dhb)-(pdy-hdy)/(dhd*dhd*dhd)-(pdy-hey)/(dhe*dhe*dhe)-(pdy-hgy)/(dhg*dhg*dhg));
     vdz=vdz+1.0e-28*ai*(-(n*pdz)/(kd*kd*kd)+(pdz-paz)/(kad*kad*kad)+(pdz-pbz)/(kbd*kbd*kbd)+(pdz-pcz)/(kcd*kcd*kcd)+(pdz-pez)/(kde*kde*kde)+(pdz-pfz)/(kdf*kdf*kdf)+(pdz-pgz)/(kdg*kdg*kdg)+(pdz-phz)/(kdh*kdh*kdh)-(pdz-hbz)/(dhb*dhb*dhb)-(pdz-hdz)/(dhd*dhd*dhd)-(pdz-hez)/(dhe*dhe*dhe)-(pdz-hgz)/(dhg*dhg*dhg));                                
    
   vex=vex+1.0e-28*ai*(-(n*pex)/(ke*ke*ke)+(pex-pax)/(kae*kae*kae)+(pex-pbx)/(kbe*kbe*kbe)+(pex-pcx)/(kce*kce*kce)+(pex-pdx)/(kde*kde*kde)+(pex-pfx)/(kef*kef*kef)+(pex-pgx)/(keg*keg*keg)+(pex-phx)/(keh*keh*keh)-(pex-hbx)/(ehb*ehb*ehb)-(pex-hdx)/(ehd*ehd*ehd)-(pex-hex)/(ehe*ehe*ehe)-(pex-hgx)/(ehg*ehg*ehg));  
     vey=vey+1.0e-28*ai*(-(n*pey)/(ke*ke*ke)+(pey-pay)/(kae*kae*kae)+(pey-pby)/(kbe*kbe*kbe)+(pey-pcy)/(kce*kce*kce)+(pey-pdy)/(kde*kde*kde)+(pey-pfy)/(kef*kef*kef)+(pey-pgy)/(keg*keg*keg)+(pey-phy)/(keh*keh*keh)-(pey-hby)/(ehb*ehb*ehb)-(pey-hdy)/(ehd*ehd*ehd)-(pey-hey)/(ehe*ehe*ehe)-(pey-hgy)/(ehg*ehg*ehg));
     vez=vez+1.0e-28*ai*(-(n*pez)/(ke*ke*ke)+(pez-paz)/(kae*kae*kae)+(pez-pbz)/(kbe*kbe*kbe)+(pez-pcz)/(kce*kce*kce)+(pez-pdz)/(kde*kde*kde)+(pez-pfz)/(kef*kef*kef)+(pez-pgz)/(keg*keg*keg)+(pez-phz)/(keh*keh*keh)-(pez-hbz)/(ehb*ehb*ehb)-(pez-hdz)/(ehd*ehd*ehd)-(pez-hez)/(ehe*ehe*ehe)-(pez-hgz)/(ehg*ehg*ehg));

   vfx=vfx+1.0e-28*ai*(-(n*pfx)/(kf*kf*kf)+(pfx-pax)/(kaf*kaf*kaf)+(pfx-pbx)/(kbf*kbf*kbf)+(pfx-pcx)/(kcf*kcf*kcf)+(pfx-pdx)/(kdf*kdf*kdf)+(pfx-pex)/(kef*kef*kef)+(pfx-pgx)/(kfg*kfg*kfg)+(pfx-phx)/(kfh*kfh*kfh)-(pfx-hbx)/(fhb*fhb*fhb)-(pfx-hdx)/(fhd*fhd*fhd)-(pfx-hex)/(fhe*fhe*fhe)-(pfx-hgx)/(fhg*fhg*fhg));  
     vfy=vfy+1.0e-28*ai*(-(n*pfy)/(kf*kf*kf)+(pfy-pay)/(kaf*kaf*kaf)+(pfy-pby)/(kbf*kbf*kbf)+(pfy-pcy)/(kcf*kcf*kcf)+(pfy-pdy)/(kdf*kdf*kdf)+(pfy-pey)/(kef*kef*kef)+(pfy-pgy)/(kfg*kfg*kfg)+(pfy-phy)/(kfh*kfh*kfh)-(pfy-hby)/(fhb*fhb*fhb)-(pfy-hdy)/(fhd*fhd*fhd)-(pfy-hey)/(fhe*fhe*fhe)-(pfy-hgy)/(fhg*fhg*fhg));
     vfz=vfz+1.0e-28*ai*(-(n*pfz)/(kf*kf*kf)+(pfz-paz)/(kaf*kaf*kaf)+(pfz-pbz)/(kbf*kbf*kbf)+(pfz-pcz)/(kcf*kcf*kcf)+(pfz-pdz)/(kdf*kdf*kdf)+(pfz-pez)/(kef*kef*kef)+(pfz-pgz)/(kfg*kfg*kfg)+(pfz-phz)/(kfh*kfh*kfh)-(pfz-hbz)/(fhb*fhb*fhb)-(pfz-hdz)/(fhd*fhd*fhd)-(pfz-hez)/(fhe*fhe*fhe)-(pfz-hgz)/(fhg*fhg*fhg));

   vgx=vgx+1.0e-28*ai*(-(n*pgx)/(kg*kg*kg)+(pgx-pax)/(kag*kag*kag)+(pgx-pbx)/(kbg*kbg*kbg)+(pgx-pcx)/(kcg*kcg*kcg)+(pgx-pdx)/(kdg*kdg*kdg)+(pgx-pex)/(keg*keg*keg)+(pgx-pfx)/(kfg*kfg*kfg)+(pgx-phx)/(kgh*kgh*kgh)-(pgx-hbx)/(ghb*ghb*ghb)-(pgx-hdx)/(ghd*ghd*ghd)-(pgx-hex)/(ghe*ghe*ghe)-(pgx-hgx)/(ghg*ghg*ghg));  
     vgy=vgy+1.0e-28*ai*(-(n*pgy)/(kg*kg*kg)+(pgy-pay)/(kag*kag*kag)+(pgy-pby)/(kbg*kbg*kbg)+(pgy-pcy)/(kcg*kcg*kcg)+(pgy-pdy)/(kdg*kdg*kdg)+(pgy-pey)/(keg*keg*keg)+(pgy-pfy)/(kfg*kfg*kfg)+(pgy-phy)/(kgh*kgh*kgh)-(pgy-hby)/(ghb*ghb*ghb)-(pgy-hdy)/(ghd*ghd*ghd)-(pgy-hey)/(ghe*ghe*ghe)-(pgy-hgy)/(ghg*ghg*ghg));
     vgz=vgz+1.0e-28*ai*(-(n*pgz)/(kg*kg*kg)+(pgz-paz)/(kag*kag*kag)+(pgz-pbz)/(kbg*kbg*kbg)+(pgz-pcz)/(kcg*kcg*kcg)+(pgz-pdz)/(kdg*kdg*kdg)+(pgz-pez)/(keg*keg*keg)+(pgz-pfz)/(kfg*kfg*kfg)+(pgz-phz)/(kgh*kgh*kgh)-(pgz-hbz)/(ghb*ghb*ghb)-(pgz-hdz)/(ghd*ghd*ghd)-(pgz-hez)/(ghe*ghe*ghe)-(pgz-hgz)/(ghg*ghg*ghg));
  
   vhx=vhx+1.0e-28*ai*(-(n*phx)/(kh*kh*kh)+(phx-pax)/(kah*kah*kah)+(phx-pbx)/(kbh*kbh*kbh)+(phx-pcx)/(kch*kch*kch)+(phx-pdx)/(kdh*kdh*kdh)+(phx-pex)/(keh*keh*keh)+(phx-pfx)/(kfh*kfh*kfh)+(phx-pgx)/(kgh*kgh*kgh)-(phx-hbx)/(hhb*hhb*hhb)-(phx-hdx)/(hhd*hhd*hhd)-(phx-hex)/(hhe*hhe*hhe)-(phx-hgx)/(hhg*hhg*hhg));  
     vhy=vhy+1.0e-28*ai*(-(n*phy)/(kh*kh*kh)+(phy-pay)/(kah*kah*kah)+(phy-pby)/(kbh*kbh*kbh)+(phy-pcy)/(kch*kch*kch)+(phy-pdy)/(kdh*kdh*kdh)+(phy-pey)/(keh*keh*keh)+(phy-pfy)/(kfh*kfh*kfh)+(phy-pgy)/(kgh*kgh*kgh)-(phy-hby)/(hhb*hhb*hhb)-(phy-hdy)/(hhd*hhd*hhd)-(phy-hey)/(hhe*hhe*hhe)-(phy-hgy)/(hhg*hhg*hhg));
     vhz=vhz+1.0e-28*ai*(-(n*phz)/(kh*kh*kh)+(phz-paz)/(kah*kah*kah)+(phz-pbz)/(kbh*kbh*kbh)+(phz-pcz)/(kch*kch*kch)+(phz-pdz)/(kdh*kdh*kdh)+(phz-pez)/(keh*keh*keh)+(phz-pfz)/(kfh*kfh*kfh)+(phz-pgz)/(kgh*kgh*kgh)-(phz-hbz)/(hhb*hhb*hhb)-(phz-hdz)/(hhd*hhd*hhd)-(phz-hez)/(hhe*hhe*hhe)-(phz-hgz)/(hhg*hhg*hhg));

    pax=ax;pay=ay;paz=az; pbx=bx;pby=by;pbz=bz; pcx=cx;pcy=cy;pcz=cz; pdx=dx;pdy=dy;pdz=dz;
    pex=ex;pey=ey;pez=ez; pfx=fx;pfy=fy;pfz=fz; pgx=gx;pgy=gy;pgz=gz; phx=hx;phy=hy;phz=hz;
  
   } while (hhrr<0);             // de Broglie's wave, wa and wb become 2.0 and 1.0 respectively           

                                // show results
  
  System.out.printf("biggest distance:%.2f ", bigdis);
  System.out.printf("shortest distance:%.2f\n", shortdis); 
   }  
   }}