healpix_base2.cc

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Copyright (C) 2003, 2004, 2005, 2006 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include "healpix_base2.h"
00033 #include "cxxutils.h"
00034 #include "pointing.h"
00035 #include "arr.h"
00036 #include "geom_utils.h"
00037 
00038 using namespace std;
00039 
00040 short Healpix_Base2::ctab[];
00041 short Healpix_Base2::utab[];
00042 
00043 Healpix_Base2::Tablefiller::Tablefiller()
00044   {
00045   for (int m=0; m<0x100; ++m)
00046     {
00047     ctab[m] =
00048          (m&0x1 )       | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
00049       | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4);
00050     utab[m] =
00051          (m&0x1 )       | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
00052       | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7);
00053     }
00054   }
00055 
00056 Healpix_Base2::Tablefiller Healpix_Base2::Filler;
00057 
00058 const int Healpix_Base2::jrll[] = { 2,2,2,2,3,3,3,3,4,4,4,4 };
00059 const int Healpix_Base2::jpll[] = { 1,3,5,7,0,2,4,6,1,3,5,7 };
00060 
00061 int64 Healpix_Base2::npix2nside (int64 npix)
00062   {
00063   int64 res=isqrt(npix/12);
00064   planck_assert (npix==res*res*12, "npix2nside: invalid argument");
00065   return res;
00066   }
00067 
00068 void Healpix_Base2::nest2xyf (int64 pix, int &ix, int &iy, int &face_num)
00069   const
00070   {
00071   face_num = pix>>(2*order_);
00072   pix &= (npface_-1);
00073   int32 raw = ((pix&0x555500000000ull)>>16) 
00074              | ((pix&0x5555000000000000ull)>>31)
00075              | (pix&0x5555)
00076              | ((pix&0x55550000)>>15);
00077   ix =  ctab[raw&0xff]
00078      | (ctab[(raw>>8)&0xff]<<4)
00079      | (ctab[(raw>>16)&0xff]<<16)
00080      | (ctab[(raw>>24)&0xff]<<20);
00081   pix >>= 1;
00082   raw = ((pix&0x555500000000ull)>>16) 
00083              | ((pix&0x5555000000000000ull)>>31)
00084              | (pix&0x5555)
00085              | ((pix&0x55550000)>>15);
00086   iy =  ctab[raw&0xff]
00087      | (ctab[(raw>>8)&0xff]<<4)
00088      | (ctab[(raw>>16)&0xff]<<16)
00089      | (ctab[(raw>>24)&0xff]<<20);
00090   }
00091 
00092 int64 Healpix_Base2::xyf2nest (int ix, int iy, int face_num) const
00093   {
00094   return (int64(face_num)<<(2*order_)) +
00095     (   (int64(utab[ ix     &0xff]))
00096       | (int64(utab[(ix>> 8)&0xff])<<16)
00097       | (int64(utab[(ix>>16)&0xff])<<32)
00098       | (int64(utab[(ix>>24)&0xff])<<48)
00099       | (int64(utab[ iy     &0xff])<<1)
00100       | (int64(utab[(iy>> 8)&0xff])<<17)
00101       | (int64(utab[(iy>>16)&0xff])<<33)
00102       | (int64(utab[(iy>>24)&0xff])<<49) ); 
00103   }
00104 
00105 void Healpix_Base2::ring2xyf (int64 pix, int &ix, int &iy, int &face_num)
00106   const
00107   {
00108   int64 iring, iphi, kshift, nr;
00109 
00110   int64 nl2 = 2*nside_;
00111 
00112   if (pix<ncap_) // North Polar cap
00113     {
00114     iring = int64(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00115     iphi  = (pix+1) - 2*iring*(iring-1);
00116     kshift = 0;
00117     nr = iring;
00118     face_num=0;
00119     int64 tmp = iphi-1;
00120     if (tmp>=(2*iring))
00121       {
00122       face_num=2;
00123       tmp-=2*iring;
00124       }
00125     if (tmp>=iring) ++face_num;
00126     }
00127   else if (pix<(npix_-ncap_)) // Equatorial region
00128     {
00129     int64 ip = pix - ncap_;
00130     if (order_>=0)
00131       {
00132       iring = (ip>>(order_+2)) + nside_; // counted from North pole
00133       iphi  = (ip&(4*nside_-1)) + 1;
00134       }
00135     else
00136       {
00137       iring = (ip/(4*nside_)) + nside_; // counted from North pole
00138       iphi  = (ip%(4*nside_)) + 1;
00139       }
00140     kshift = (iring+nside_)&1;
00141     nr = nside_;
00142     int64 ire = iring-nside_+1;
00143     int64 irm = nl2+2-ire;
00144     int64 ifm, ifp;
00145     if (order_>=0)
00146       {
00147       ifm = (iphi - ire/2 + nside_ -1) >> order_;
00148       ifp = (iphi - irm/2 + nside_ -1) >> order_;
00149       }
00150     else
00151       {
00152       ifm = (iphi - ire/2 + nside_ -1) / nside_;
00153       ifp = (iphi - irm/2 + nside_ -1) / nside_;
00154       }
00155     if (ifp == ifm) // faces 4 to 7
00156       face_num = (ifp==4) ? 4 : ifp+4;
00157     else if (ifp<ifm) // (half-)faces 0 to 3
00158       face_num = ifp;
00159     else // (half-)faces 8 to 11
00160       face_num = ifm + 8;
00161     }
00162   else // South Polar cap
00163     {
00164     int64 ip = npix_ - pix;
00165     iring = int64(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00166     iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
00167     kshift = 0;
00168     nr = iring;
00169     iring = 2*nl2-iring;
00170     face_num=8;
00171     int64 tmp = iphi-1;
00172     if (tmp>=(2*nr))
00173       {
00174       face_num=10;
00175       tmp-=2*nr;
00176       }
00177     if (tmp>=nr) ++face_num;
00178     }
00179 
00180   int64 irt = iring - (jrll[face_num]*nside_) + 1;
00181   int64 ipt = 2*iphi- jpll[face_num]*nr - kshift -1;
00182   if (ipt>=nl2) ipt-=8*nside_;
00183 
00184   ix =  (ipt-irt) >>1;
00185   iy =(-(ipt+irt))>>1;
00186   }
00187 
00188 int64 Healpix_Base2::xyf2ring (int ix, int iy, int face_num) const
00189   {
00190   int64 nl4 = 4*nside_;
00191   int64 jr = (jrll[face_num]*nside_) - ix - iy  - 1;
00192 
00193   int64 nr, kshift, n_before;
00194   if (jr<nside_)
00195     {
00196     nr = jr;
00197     n_before = 2*nr*(nr-1);
00198     kshift = 0;
00199     }
00200   else if (jr > 3*nside_)
00201     {
00202     nr = nl4-jr;
00203     n_before = npix_ - 2*(nr+1)*nr;
00204     kshift = 0;
00205     }
00206   else
00207     {
00208     nr = nside_;
00209     n_before = ncap_ + (jr-nside_)*nl4;
00210     kshift = (jr-nside_)&1;
00211     }
00212 
00213   int64 jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
00214   if (jp>nl4)
00215     jp-=nl4;
00216   else
00217     if (jp<1) jp+=nl4;
00218 
00219   return n_before + jp - 1;
00220   }
00221 
00222 double Healpix_Base2::ring2z (int64 ring) const
00223   {
00224   if (ring<nside_)
00225     return 1 - ring*ring*fact2_;
00226   if (ring <=3*nside_)
00227     return (2*nside_-ring)*fact1_;
00228   ring=4*nside_ - ring;
00229   return ring*ring*fact2_ - 1;
00230   }
00231 
00232 int64 Healpix_Base2::pix2ring (int64 pix) const
00233   {
00234   if (scheme_==RING)
00235     {
00236     if (pix<ncap_) // North Polar cap
00237       return int64(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00238     else if (pix<(npix_-ncap_)) // Equatorial region
00239       {
00240       int64 ip  = pix - ncap_;
00241       return ip/(4*nside_) + nside_; // counted from North pole
00242       }
00243     else // South Polar cap
00244       {
00245       int64 ip = npix_ - pix;
00246       return 4*nside_ - int64(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00247       }
00248     }
00249   else
00250     {
00251     int face_num, ix, iy;
00252     nest2xyf(pix,ix,iy,face_num);
00253     return (int64(jrll[face_num])<<order_) - ix - iy - 1;
00254     }
00255   }
00256 
00257 int64 Healpix_Base2::nest2ring (int64 pix) const
00258   {
00259   planck_assert(order_>=0, "nest2ring: need hierarchical map");
00260   int ix, iy, face_num;
00261   nest2xyf (pix, ix, iy, face_num);
00262   return xyf2ring (ix, iy, face_num);
00263   }
00264 
00265 int64 Healpix_Base2::ring2nest (int64 pix) const
00266   {
00267   planck_assert(order_>=0, "ring2nest: need hierarchical map");
00268   int ix, iy, face_num;
00269   ring2xyf (pix, ix, iy, face_num);
00270   return xyf2nest (ix, iy, face_num);
00271   }
00272 
00273 int64 Healpix_Base2::nest2peano (int64 pix) const
00274   {
00275   static const unsigned char subpix[8][4] = {
00276     { 0, 1, 3, 2 }, { 3, 0, 2, 1 }, { 2, 3, 1, 0 }, { 1, 2, 0, 3 },
00277     { 0, 3, 1, 2 }, { 1, 0, 2, 3 }, { 2, 1, 3, 0 }, { 3, 2, 0, 1 } };
00278   const unsigned char subpath[8][4] = {
00279     { 4, 0, 6, 0 }, { 7, 5, 1, 1 }, { 2, 4, 2, 6 }, { 3, 3, 7, 5 },
00280     { 0, 2, 4, 4 }, { 5, 1, 5, 3 }, { 6, 6, 0, 2 }, { 1, 7, 3, 7 } };
00281   static const unsigned char face2path[12] = {
00282     2, 5, 2, 5, 3, 6, 3, 6, 2, 3, 2, 3 };
00283   static const unsigned char face2peanoface[12] = {
00284     0, 5, 6, 11, 10, 1, 4, 7, 2, 3, 8, 9 };
00285 
00286   int face = pix>>(2*order_);
00287   unsigned char path = face2path[face];
00288   int64 result = 0;
00289 
00290   for (int shift=2*order_-2; shift>=0; shift-=2)
00291     {
00292     unsigned char spix = (pix>>shift) & 0x3;
00293     result <<= 2;
00294     result |= subpix[path][spix];
00295     path=subpath[path][spix];
00296     }
00297 
00298   return result + (int64(face2peanoface[face])<<(2*order_));
00299   }
00300 
00301 int64 Healpix_Base2::peano2nest (int64 pix) const
00302   {
00303   static const unsigned char subpix[8][4] = {
00304     { 0, 1, 3, 2 }, { 1, 3, 2, 0 }, { 3, 2, 0, 1 }, { 2, 0, 1, 3 },
00305     { 0, 2, 3, 1 }, { 1, 0, 2, 3 }, { 3, 1, 0, 2 }, { 2, 3, 1, 0 } };
00306   static const unsigned char subpath[8][4] = {
00307     { 4, 0, 0, 6 }, { 5, 1, 1, 7 }, { 6, 2, 2, 4 }, { 7, 3, 3, 5 },
00308     { 0, 4, 4, 2 }, { 1, 5, 5, 3 }, { 2, 6, 6, 0 }, { 3, 7, 7, 1 } };
00309   static const unsigned char face2path[12] = {
00310     2, 6, 2, 3, 3, 5, 2, 6, 2, 3, 3, 5 };
00311   static const unsigned char peanoface2face[12] = {
00312     0, 5, 8, 9, 6, 1, 2, 7, 10, 11, 4, 3 };
00313 
00314   int face = pix>>(2*order_);
00315   unsigned char path = face2path[face];
00316   int64 result = 0;
00317 
00318   for (int shift=2*order_-2; shift>=0; shift-=2)
00319     {
00320     unsigned char spix = (pix>>shift) & 0x3;
00321     result <<= 2;
00322     result |= subpix[path][spix];
00323     path=subpath[path][spix];
00324     }
00325 
00326   return result + (int64(peanoface2face[face])<<(2*order_));
00327   }
00328 
00329 int64 Healpix_Base2::ang2pix_z_phi (double z, double phi) const
00330   {
00331   double za = abs(z);
00332   double tt = fmodulo(phi,twopi) * inv_halfpi; // in [0,4)
00333 
00334   if (scheme_==RING)
00335     {
00336     if (za<=twothird) // Equatorial region
00337       {
00338       double temp1 = nside_*(0.5+tt);
00339       double temp2 = nside_*z*0.75;
00340       int64 jp = int64(temp1-temp2); // index of  ascending edge line
00341       int64 jm = int64(temp1+temp2); // index of descending edge line
00342 
00343       // ring number counted from z=2/3
00344       int64 ir = nside_ + 1 + jp - jm; // in {1,2n+1}
00345       int kshift = 1-(ir&1); // kshift=1 if ir even, 0 otherwise
00346 
00347       int64 ip = (jp+jm-nside_+kshift+1)/2; // in {0,4n-1}
00348       ip = imodulo<int64>(ip,4*nside_);
00349 
00350       return ncap_ + (ir-1)*4*nside_ + ip;
00351       }
00352     else  // North & South polar caps
00353       {
00354       double tp = tt-int(tt);
00355       double tmp = nside_*sqrt(3*(1-za));
00356 
00357       int64 jp = int64(tp*tmp); // increasing edge line index
00358       int64 jm = int64((1.0-tp)*tmp); // decreasing edge line index
00359 
00360       int64 ir = jp+jm+1; // ring number counted from the closest pole
00361       int64 ip = int64(tt*ir); // in {0,4*ir-1}
00362       ip = imodulo<int64>(ip,4*ir);
00363 
00364       if (z>0)
00365         return 2*ir*(ir-1) + ip;
00366       else
00367         return npix_ - 2*ir*(ir+1) + ip;
00368       }
00369     }
00370   else // scheme_ == NEST
00371     {
00372     int face_num, ix, iy;
00373 
00374     if (za<=twothird) // Equatorial region
00375       {
00376       double temp1 = nside_*(0.5+tt);
00377       double temp2 = nside_*(z*0.75);
00378       int64 jp = int64(temp1-temp2); // index of  ascending edge line
00379       int64 jm = int64(temp1+temp2); // index of descending edge line
00380       int64 ifp = jp >> order_;  // in {0,4}
00381       int64 ifm = jm >> order_;
00382       if (ifp == ifm)           // faces 4 to 7
00383         face_num = (ifp==4) ? 4: ifp+4;
00384       else if (ifp < ifm)       // (half-)faces 0 to 3
00385         face_num = ifp;
00386       else                      // (half-)faces 8 to 11
00387         face_num = ifm + 8;
00388 
00389       ix = jm & (nside_-1);
00390       iy = nside_ - (jp & (nside_-1)) - 1;
00391       }
00392     else // polar region, za > 2/3
00393       {
00394       int ntt = int(tt);
00395       if (ntt>=4) ntt=3;
00396       double tp = tt-ntt;
00397       double tmp = nside_*sqrt(3*(1-za));
00398 
00399       int64 jp = int64(tp*tmp); // increasing edge line index
00400       int64 jm = int64((1.0-tp)*tmp); // decreasing edge line index
00401       if (jp>=nside_) jp = nside_-1; // for points too close to the boundary
00402       if (jm>=nside_) jm = nside_-1;
00403       if (z >= 0)
00404         {
00405         face_num = ntt;  // in {0,3}
00406         ix = nside_ - jm - 1;
00407         iy = nside_ - jp - 1;
00408         }
00409       else
00410         {
00411         face_num = ntt + 8; // in {8,11}
00412         ix =  jp;
00413         iy =  jm;
00414         }
00415       }
00416 
00417     return xyf2nest(ix,iy,face_num);
00418     }
00419   }
00420 
00421 void Healpix_Base2::pix2ang_z_phi (int64 pix, double &z, double &phi) const
00422   {
00423   if (scheme_==RING)
00424     {
00425     if (pix<ncap_) // North Polar cap
00426       {
00427       int64 iring = int64(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00428       int64 iphi  = (pix+1) - 2*iring*(iring-1);
00429 
00430       z = 1.0 - (iring*iring)*fact2_;
00431       phi = (iphi-0.5) * halfpi/iring;
00432       }
00433     else if (pix<(npix_-ncap_)) // Equatorial region
00434       {
00435       int64 ip  = pix - ncap_;
00436       int64 iring = ip/(4*nside_) + nside_; // counted from North pole
00437       int64 iphi  = ip%(4*nside_) + 1;
00438       // 1 if iring+nside is odd, 1/2 otherwise
00439       double fodd = ((iring+nside_)&1) ? 1 : 0.5;
00440 
00441       int64 nl2 = 2*nside_;
00442       z = (nl2-iring)*fact1_;
00443       phi = (iphi-fodd) * pi/nl2;
00444       }
00445     else // South Polar cap
00446       {
00447       int64 ip = npix_ - pix;
00448       int64 iring = int64(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00449       int64 iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
00450 
00451       z = -1.0 + (iring*iring)*fact2_;
00452       phi = (iphi-0.5) * halfpi/iring;
00453       }
00454     }
00455   else
00456     {
00457     int64 nl4 = nside_*4;
00458 
00459     int face_num, ix, iy;
00460     nest2xyf(pix,ix,iy,face_num);
00461 
00462     int64 jr = (int64(jrll[face_num])<<order_) - ix - iy - 1;
00463 
00464     int64 nr;
00465     int kshift;
00466     if (jr<nside_)
00467       {
00468       nr = jr;
00469       z = 1 - nr*nr*fact2_;
00470       kshift = 0;
00471       }
00472     else if (jr > 3*nside_)
00473       {
00474       nr = nl4-jr;
00475       z = nr*nr*fact2_ - 1;
00476       kshift = 0;
00477       }
00478     else
00479       {
00480       nr = nside_;
00481       z = (2*nside_-jr)*fact1_;
00482       kshift = (jr-nside_)&1;
00483       }
00484 
00485     int64 jp = (jpll[face_num]*nr + ix -iy + 1 + kshift) / 2;
00486     if (jp>nl4) jp-=nl4;
00487     if (jp<1) jp+=nl4;
00488 
00489     phi = (jp-(kshift+1)*0.5)*(halfpi/nr);
00490     }
00491   }
00492 
00493 void Healpix_Base2::neighbors (int64 pix, fix_arr<int64,8> &result) const
00494   {
00495   static const int xoffset[] = { -1,-1, 0, 1, 1, 1, 0,-1 };
00496   static const int yoffset[] = {  0, 1, 1, 1, 0,-1,-1,-1 };
00497   static const int facearray[][12] =
00498         { {  8, 9,10,11,-1,-1,-1,-1,10,11, 8, 9 },   // S
00499           {  5, 6, 7, 4, 8, 9,10,11, 9,10,11, 8 },   // SE
00500           { -1,-1,-1,-1, 5, 6, 7, 4,-1,-1,-1,-1 },   // E
00501           {  4, 5, 6, 7,11, 8, 9,10,11, 8, 9,10 },   // SW
00502           {  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11 },   // center
00503           {  1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 },   // NE
00504           { -1,-1,-1,-1, 7, 4, 5, 6,-1,-1,-1,-1 },   // W
00505           {  3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 },   // NW
00506           {  2, 3, 0, 1,-1,-1,-1,-1, 0, 1, 2, 3 } }; // N
00507   static const int swaparray[][12] =
00508         { {  0,0,0,0,0,0,0,0,3,3,3,3 },   // S
00509           {  0,0,0,0,0,0,0,0,6,6,6,6 },   // SE
00510           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // E
00511           {  0,0,0,0,0,0,0,0,5,5,5,5 },   // SW
00512           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // center
00513           {  5,5,5,5,0,0,0,0,0,0,0,0 },   // NE
00514           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // W
00515           {  6,6,6,6,0,0,0,0,0,0,0,0 },   // NW
00516           {  3,3,3,3,0,0,0,0,0,0,0,0 } }; // N
00517 
00518   int ix, iy, face_num;
00519   (scheme_==RING) ?
00520     ring2xyf(pix,ix,iy,face_num) : nest2xyf(pix,ix,iy,face_num);
00521 
00522   const int64 nsm1 = nside_-1;
00523   if ((ix>0)&&(ix<nsm1)&&(iy>0)&&(iy<nsm1))
00524     {
00525     if (scheme_==RING)
00526       for (int m=0; m<8; ++m)
00527         result[m] = xyf2ring(ix+xoffset[m],iy+yoffset[m],face_num);
00528     else
00529       for (int m=0; m<8; ++m)
00530         result[m] = xyf2nest(ix+xoffset[m],iy+yoffset[m],face_num);
00531     }
00532   else
00533     {
00534     for (int i=0; i<8; ++i)
00535       {
00536       int x=ix+xoffset[i];
00537       int y=iy+yoffset[i];
00538       int nbnum=4;
00539       if (x<0)
00540         { x+=nside_; nbnum-=1; }
00541       else if (x>=nside_)
00542         { x-=nside_; nbnum+=1; }
00543       if (y<0)
00544         { y+=nside_; nbnum-=3; }
00545       else if (y>=nside_)
00546         { y-=nside_; nbnum+=3; }
00547 
00548       int f = facearray[nbnum][face_num];
00549       if (f>=0)
00550         {
00551         if (swaparray[nbnum][face_num]&1) x=nside_-x-1;
00552         if (swaparray[nbnum][face_num]&2) y=nside_-y-1;
00553         if (swaparray[nbnum][face_num]&4) std::swap(x,y);
00554         result[i] = (scheme_==RING) ? xyf2ring(x,y,f) : xyf2nest(x,y,f);
00555         }
00556       else
00557         result[i] = -1;
00558       }
00559     }
00560   }
00561 
00562 double Healpix_Base2::max_pixrad() const
00563   {
00564   vec3 va,vb;
00565   va.set_z_phi (2./3., pi/(4*nside_));
00566   double t1 = 1.-1./nside_;
00567   t1*=t1;
00568   vb.set_z_phi (1-t1/3, 0);
00569   return v_angle(va,vb);
00570   }