ADMB Documentation  11.5.3152
 All Classes Files Functions Variables Typedefs Friends Defines
df33fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00007 #include <df1b2fun.h>
00008 #include "df33fun.h"
00009   df1b2variable * df3_three_variable::ind_var[3];
00010   int df3_three_variable::num_ind_var=0;
00011   int df3_three_variable::num_local_ind_var=0;
00012 
00016 df3_three_variable::df3_three_variable()
00017 {
00018   for (int i = 0; i < 20; ++i)
00019   {
00020     v[i] = 0;
00021   }
00022 }
00026 df3_three_variable::df3_three_variable(const df3_three_variable& x)
00027 {
00028   memcpy(&(v[0]),&(x.v[0]),20*sizeof(double));
00029 }
00033 df3_three_vector::df3_three_vector(const df3_three_vector& m2)
00034 {
00035   index_min=m2.index_min;
00036   index_max=m2.index_max;
00037   shape=m2.shape;
00038   if (shape)
00039   {
00040     (shape->ncopies)++;
00041   }
00042   v = m2.v;
00043 }
00044 
00049 df1b2variable& df1b2variable::operator = (const df3_three_variable& v)
00050 {
00051   df1b2variable * px=df3_three_variable::ind_var[0];
00052   df1b2variable * py=df3_three_variable::ind_var[1];
00053   df1b2variable * pv=df3_three_variable::ind_var[2];
00054   df3_three_variable::num_ind_var=0;
00055   df3_three_variable::ind_var[0]=0;
00056   df3_three_variable::ind_var[1]=0;
00057   df3_three_variable::ind_var[2]=0;
00058   //df1b2variable * px=0;
00059   double  dfx= *v.get_u_x();
00060   double  dfy= *v.get_u_y();
00061   double  dfz= *v.get_u_z();
00062   double * xd=px->get_u_dot();
00063   double * yd=py->get_u_dot();
00064   double * vd=pv->get_u_dot();
00065   double * zd=get_u_dot();
00066   *get_u()=*v.get_u();
00067   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00068   {
00069     *zd++ = dfx * *xd++ + dfy * *yd++ + dfz * *vd++;
00070   }
00071 
00072  /*
00073   cout << *v.get_u()  << " ";
00074   cout << *v.get_udot()  << " ";
00075   cout << *v.get_udot2()  << " ";
00076   cout << *v.get_udot3()  << endl;
00077   */
00078   f1b2gradlist->write_pass1(px,
00079     py,
00080     pv,
00081     this,
00082     *v.get_u_x(),
00083     *v.get_u_y(),
00084     *v.get_u_z(),
00085     *v.get_u_xx(),
00086     *v.get_u_xy(),
00087     *v.get_u_xz(),
00088     *v.get_u_yy(),
00089     *v.get_u_yz(),
00090     *v.get_u_zz(),
00091     *v.get_u_xxx(),
00092     *v.get_u_xxy(),
00093     *v.get_u_xxz(),
00094     *v.get_u_xyy(),
00095     *v.get_u_xyz(),
00096     *v.get_u_xzz(),
00097     *v.get_u_yyy(),
00098     *v.get_u_yyz(),
00099     *v.get_u_yzz(),
00100     *v.get_u_zzz());
00101   return *this;
00102 }
00106 void df3_three_variable::initialize()
00107 {
00108   for (int i = 0; i < 20; ++i)
00109   {
00110     v[i] = 0.0;
00111   }
00112 }
00116 df3_three_vector::~df3_three_vector()
00117 {
00118   deallocate();
00119 }
00120 
00125  void df3_three_vector::deallocate(void)
00126  {
00127    if(shape)
00128    {
00129      if (shape->ncopies)
00130      {
00131        (shape->ncopies)--;
00132      }
00133      else
00134      {
00135        v = (df3_three_variable*) (shape->trueptr);
00136        delete [] v;
00137        v = NULL;
00138        delete shape;
00139        shape=0;
00140      }
00141    }
00142  }
00143 
00148  dvector value(const df3_three_vector& v)
00149  {
00150    int mmin=v.indexmin();
00151    int mmax=v.indexmax();
00152    dvector cv(mmin,mmax);
00153    for (int i=mmin;i<=mmax;i++)
00154    {
00155      cv(i)=value(v(i));
00156    }
00157    return cv;
00158  }
00159 
00164   void df3_three_vector::initialize(void)
00165   {
00166     int mmin=indexmin();
00167     int mmax=indexmax();
00168     for (int i=mmin;i<=mmax;i++)
00169     {
00170       (*this)(i)=0.0;
00171     }
00172   }
00176 df3_three_vector::df3_three_vector()
00177 {
00178   allocate();
00179 }
00183 df3_three_vector::df3_three_vector(int min,int max)
00184 {
00185   allocate(min, max);
00186 }
00187 
00192   void df3_three_vector::allocate(int min,int max)
00193   {
00194     index_min=min;
00195     index_max=max;
00196     v=new df3_three_variable[max-min+1];
00197     if (v==0)
00198     {
00199       cerr << "error allocating memory in df3_three_vector" << endl;
00200       ad_exit(1);
00201     }
00202     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00203     {
00204       cerr << "Error trying to allocate memory for df3_three_vector"
00205            << endl;;
00206       ad_exit(1);
00207     }
00208     v-=min;
00209   }
00213 void df3_three_vector::allocate(void)
00214 {
00215   index_min = 0;
00216   index_max = -1;
00217   v = NULL;
00218   shape = NULL;
00219 }
00220 
00225  dmatrix value(const df3_three_matrix& v)
00226  {
00227    int rmin=v.indexmin();
00228    int rmax=v.indexmax();
00229    dmatrix cm(rmin,rmax);
00230    for (int i=rmin;i<=rmax;i++)
00231    {
00232      int cmin=v(i).indexmin();
00233      int cmax=v(i).indexmax();
00234      cm(i).allocate(cmin,cmax);
00235      for (int j=cmin;j<=cmax;j++)
00236      {
00237        cm(i,j)=value(v(i,j));
00238      }
00239    }
00240    return cm;
00241  }
00242 
00247  df3_three_matrix::df3_three_matrix(const df3_three_matrix& m2)
00248  {
00249    index_min=m2.index_min;
00250    index_max=m2.index_max;
00251    shape=m2.shape;
00252    if (shape)
00253    {
00254      (shape->ncopies)++;
00255    }
00256    v = m2.v;
00257  }
00258 
00263  df3_three_matrix::~df3_three_matrix()
00264  {
00265    deallocate();
00266  }
00267 
00272  void df3_three_matrix::deallocate(void)
00273  {
00274    if (shape)
00275    {
00276      if (shape->ncopies)
00277      {
00278        (shape->ncopies)--;
00279      }
00280      else
00281      {
00282        v = (df3_three_vector*) (shape->get_pointer());
00283        delete [] v;
00284        v=0;
00285        delete shape;
00286        shape=0;
00287      }
00288    }
00289  }
00290 
00295   void df3_three_matrix::initialize(void)
00296   {
00297     int mmin=indexmin();
00298     int mmax=indexmax();
00299     for (int i=mmin;i<=mmax;i++)
00300     {
00301       (*this)(i).initialize();
00302     }
00303   }
00304 
00309   df3_three_matrix::df3_three_matrix(int rmin,int rmax,int cmin,int cmax)
00310   {
00311     index_min=rmin;
00312     index_max=rmax;
00313     v=new df3_three_vector[rmax-rmin+1];
00314     if (v==0)
00315     {
00316       cerr << "error allocating memory in df3_three_matrix" << endl;
00317       ad_exit(1);
00318     }
00319     if ( (shape=new mat_shapex(v)) == NULL)
00320     {
00321       cerr << "Error trying to allocate memory for df3_three_vector"
00322            << endl;;
00323     }
00324     v-=rmin;
00325 
00326     for (int i=rmin;i<=rmax;i++)
00327     {
00328       v[i].allocate(cmin,cmax);
00329     }
00330   }
00334 df3_three_variable& df3_three_variable::operator-=(const df3_three_variable& _v)
00335 {
00336   *get_u() -= *_v.get_u();
00337   *get_u_x() -= *_v.get_u_x();
00338   *get_u_y() -= *_v.get_u_y();
00339   *get_u_z() -= *_v.get_u_z();
00340   *get_u_xx() -= *_v.get_u_xx();
00341   *get_u_xy() -= *_v.get_u_xy();
00342   *get_u_xz() -= *_v.get_u_xz();
00343   *get_u_yy() -= *_v.get_u_yy();
00344   *get_u_yz() -= *_v.get_u_yz();
00345   *get_u_zz() -= *_v.get_u_zz();
00346   *get_u_xxx() -= *_v.get_u_xxx();
00347   *get_u_xxy() -= *_v.get_u_xxy();
00348   *get_u_xxz() -= *_v.get_u_xxz();
00349   *get_u_xyy() -= *_v.get_u_xyy();
00350   *get_u_xyz() -= *_v.get_u_xyz();
00351   *get_u_xzz() -= *_v.get_u_xzz();
00352   *get_u_yyy() -= *_v.get_u_yyy();
00353   *get_u_yyz() -= *_v.get_u_yyz();
00354   *get_u_yzz() -= *_v.get_u_yzz();
00355   *get_u_zzz() -= *_v.get_u_zzz();
00356   return *this;
00357 }
00361 df3_three_variable& df3_three_variable::operator-=(double _v)
00362 {
00363   *get_u() -= _v;
00364   return *this;
00365 }
00366 
00371 df3_three_variable operator-(const df3_three_variable& v)
00372 {
00373   df3_three_variable z;
00374 
00375   *z.get_u() = -(*v.get_u());
00376   *z.get_u_x() = -(*v.get_u_x());
00377   *z.get_u_y() = -(*v.get_u_y());
00378   *z.get_u_z() = -(*v.get_u_z());
00379   *z.get_u_xx() = -(*v.get_u_xx());
00380   *z.get_u_xy() = -(*v.get_u_xy());
00381   *z.get_u_xz() = -(*v.get_u_xz());
00382   *z.get_u_yy() = -(*v.get_u_yy());
00383   *z.get_u_yz() = -(*v.get_u_yz());
00384   *z.get_u_zz() = -(*v.get_u_zz());
00385   *z.get_u_xxx() = -(*v.get_u_xxx());
00386   *z.get_u_xxy() = -(*v.get_u_xxy());
00387   *z.get_u_xxz() = -(*v.get_u_xxz());
00388   *z.get_u_xyy() = -(*v.get_u_xyy());
00389   *z.get_u_xyz() = -(*v.get_u_xyz());
00390   *z.get_u_xzz() = -(*v.get_u_xzz());
00391   *z.get_u_yyy() = -(*v.get_u_yyy());
00392   *z.get_u_yyz() = -(*v.get_u_yyz());
00393   *z.get_u_yzz() = -(*v.get_u_yzz());
00394   *z.get_u_zzz() = -(*v.get_u_zzz());
00395 
00396   return z;
00397 }
00401 df3_three_variable& df3_three_variable::operator+=(const df3_three_variable& _v)
00402 {
00403   *get_u() += *_v.get_u();
00404   *get_u_x() += *_v.get_u_x();
00405   *get_u_y() += *_v.get_u_y();
00406   *get_u_z() += *_v.get_u_z();
00407   *get_u_xx() += *_v.get_u_xx();
00408   *get_u_xy() += *_v.get_u_xy();
00409   *get_u_xz() += *_v.get_u_xz();
00410   *get_u_yy() += *_v.get_u_yy();
00411   *get_u_yz() += *_v.get_u_yz();
00412   *get_u_zz() += *_v.get_u_zz();
00413   *get_u_xxx() += *_v.get_u_xxx();
00414   *get_u_xxy() += *_v.get_u_xxy();
00415   *get_u_xxz() += *_v.get_u_xxz();
00416   *get_u_xyy() += *_v.get_u_xyy();
00417   *get_u_xyz() += *_v.get_u_xyz();
00418   *get_u_xzz() += *_v.get_u_xzz();
00419   *get_u_yyy() += *_v.get_u_yyy();
00420   *get_u_yyz() += *_v.get_u_yyz();
00421   *get_u_yzz() += *_v.get_u_yzz();
00422   *get_u_zzz() += *_v.get_u_zzz();
00423   return *this;
00424 }
00428 df3_three_variable& df3_three_variable::operator*=(double _v)
00429 {
00430   *get_u() *= _v;
00431   *get_u_x() *= _v;
00432   *get_u_y() *= _v;
00433   *get_u_z() *= _v;
00434   *get_u_xx() *= _v;
00435   *get_u_xy() *= _v;
00436   *get_u_xz() *= _v;
00437   *get_u_yy() *= _v;
00438   *get_u_yz() *= _v;
00439   *get_u_zz() *= _v;
00440   *get_u_xxx() *= _v;
00441   *get_u_xxy() *= _v;
00442   *get_u_xxz() *= _v;
00443   *get_u_xyy() *= _v;
00444   *get_u_xyz() *= _v;
00445   *get_u_xzz() *= _v;
00446   *get_u_yyy() *= _v;
00447   *get_u_yyz() *= _v;
00448   *get_u_yzz() *= _v;
00449   *get_u_zzz() *= _v;
00450   return *this;
00451 }
00455 df3_three_variable& df3_three_variable::operator*=(const df3_three_variable& v)
00456 {
00457   df3_three_variable x = *this * v;
00458   *this = x;
00459   return *this;
00460 }
00464 df3_three_variable& df3_three_variable::operator/=(const df3_three_variable& v)
00465 {
00466   df3_three_variable x = *this / v;
00467   *this = x;
00468   return *this;
00469 }
00473 df3_three_variable& df3_three_variable::operator+=(double _v)
00474 {
00475   *get_u() += _v;
00476   return *this;
00477 }
00478 
00483 void set_derivatives(df3_three_variable& z,const df3_three_variable& x,
00484   double u, double zp,double zp2,double zp3)
00485 {
00486     //*z.get_u() = u;
00487 
00488     *z.get_u_x() = zp* *x.get_u_x();
00489 
00490     *z.get_u_y() = zp* *x.get_u_y();
00491 
00492     *z.get_u_z() = zp* *x.get_u_z();
00493 
00494     *z.get_u_xx() = zp2 * square(*x.get_u_x())
00495                    + zp * *x.get_u_xx();
00496 
00497     *z.get_u_xy() = zp2 * *x.get_u_x() * *x.get_u_y()
00498                    + zp * *x.get_u_xy();
00499 
00500     *z.get_u_xz() = zp2 * *x.get_u_x() * *x.get_u_z()
00501                    + zp * *x.get_u_xz();
00502 
00503     *z.get_u_yy() = zp2 * square(*x.get_u_y())
00504                    + zp * *x.get_u_yy();
00505 
00506     *z.get_u_yz() = zp2 * *x.get_u_y()* *x.get_u_z()
00507                    + zp * *x.get_u_yz();
00508 
00509     *z.get_u_zz() = zp2 * square(*x.get_u_z())
00510                    + zp * *x.get_u_zz();
00511 
00512     *z.get_u_xxx() = zp3 * cube(*x.get_u_x())
00513                    + 3.0 * zp2 * *x.get_u_x() * *x.get_u_xx()
00514                    + zp * *x.get_u_xxx();
00515 
00516     *z.get_u_xxy() = zp3 * square(*x.get_u_x())* *x.get_u_y()
00517                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xy()
00518                    + zp2 * *x.get_u_y() * *x.get_u_xx()
00519                    + zp * *x.get_u_xxy();
00520 
00521     *z.get_u_xxz() = zp3 * square(*x.get_u_x())* *x.get_u_z()
00522                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xz()
00523                    + zp2 * *x.get_u_z() * *x.get_u_xx()
00524                    + zp * *x.get_u_xxz();
00525 
00526     *z.get_u_xyy() = zp3 * square(*x.get_u_y())* *x.get_u_x()
00527                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_xy()
00528                    + zp2 * *x.get_u_x() * *x.get_u_yy()
00529                    + zp * *x.get_u_xyy();
00530 
00531     *z.get_u_xyz() = zp3 * *x.get_u_y()* *x.get_u_x()* *x.get_u_z()
00532                    + zp2 * *x.get_u_z() * *x.get_u_xy()
00533                    + zp2 * *x.get_u_y() * *x.get_u_xz()
00534                    + zp2 * *x.get_u_x() * *x.get_u_yz()
00535                    + zp * *x.get_u_xyz();
00536 
00537     *z.get_u_xzz() = zp3 * square(*x.get_u_z())* *x.get_u_x()
00538                    + 2.0* zp2 * *x.get_u_z() * *x.get_u_xz()
00539                    + zp2 * *x.get_u_x() * *x.get_u_zz()
00540                    + zp * *x.get_u_xzz();
00541 
00542 
00543     *z.get_u_yyy() = zp3 * cube(*x.get_u_y())
00544                    + 3.0 * zp2 * *x.get_u_y() * *x.get_u_yy()
00545                    + zp * *x.get_u_yyy();
00546     *z.get_u_yyz() = zp3 * square(*x.get_u_y())* *x.get_u_z()
00547                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_yz()
00548                    + zp2 * *x.get_u_z() * *x.get_u_yy()
00549                    + zp * *x.get_u_yyz();
00550 
00551     *z.get_u_yzz() = zp3 * square(*x.get_u_z())* *x.get_u_y()
00552                    + 2.0* zp2 * *x.get_u_z() * *x.get_u_yz()
00553                    + zp2 * *x.get_u_y() * *x.get_u_zz()
00554                    + zp * *x.get_u_yzz();
00555 
00556     *z.get_u_zzz() = zp3 * cube(*x.get_u_z())
00557                    + 3.0 * zp2 * *x.get_u_z() * *x.get_u_zz()
00558                    + zp * *x.get_u_zzz();
00559 }
00560 
00565 void set_derivatives( df3_three_variable& z, const df3_three_variable& x,
00566   const df3_three_variable& y, double u,
00567   double f_u,double f_v,double f_uu,double f_uv,double f_vv,double f_uuu,
00568   double f_uuv,double f_uvv,double f_vvv)
00569 {
00570     *z.get_u() = u;
00571 
00572     *z.get_u_x() = f_u* *x.get_u_x()
00573                  + f_v* *y.get_u_x();
00574 
00575     *z.get_u_y() = f_u* *x.get_u_y()
00576                  + f_v* *y.get_u_y();
00577 
00578     *z.get_u_z() = f_u* *x.get_u_z()
00579                  + f_v* *y.get_u_z();
00580 
00581     *z.get_u_xx() = f_uu * square(*x.get_u_x())
00582                   + f_u  * *x.get_u_xx()
00583                   + f_vv * square(*y.get_u_x())
00584                   + f_v  * *y.get_u_xx()
00585             + 2.0 * f_uv * *x.get_u_x() * *y.get_u_x();
00586 
00587     *z.get_u_xy() = f_uu * *x.get_u_x() * *x.get_u_y()
00588                   + f_u  * *x.get_u_xy()
00589                   + f_vv *  *y.get_u_x() * *y.get_u_y()
00590                   + f_v  * *y.get_u_xy()
00591                   + f_uv * (*x.get_u_x() * *y.get_u_y()
00592                          +  *x.get_u_y() * *y.get_u_x());
00593 
00594     *z.get_u_xz() = f_uu * *x.get_u_x() * *x.get_u_z()
00595                   + f_u  * *x.get_u_xz()
00596                   + f_vv *  *y.get_u_x() * *y.get_u_z()
00597                   + f_v  * *y.get_u_xz()
00598                   + f_uv * (*x.get_u_x() * *y.get_u_z()
00599                          +  *x.get_u_z() * *y.get_u_x());
00600 
00601     *z.get_u_yy() = f_uu * square(*x.get_u_y())
00602                   + f_u  * *x.get_u_yy()
00603                   + f_vv * square(*y.get_u_y())
00604                   + f_v  * *y.get_u_yy()
00605             + 2.0 * f_uv * *x.get_u_y() * *y.get_u_y();
00606 
00607     *z.get_u_yz() = f_uu * *x.get_u_y() * *x.get_u_z()
00608                   + f_u  * *x.get_u_yz()
00609                   + f_vv *  *y.get_u_y() * *y.get_u_z()
00610                   + f_v  * *y.get_u_yz()
00611                   + f_uv * (*x.get_u_y() * *y.get_u_z()
00612                          +  *x.get_u_z() * *y.get_u_y());
00613 
00614     *z.get_u_zz() = f_uu * *x.get_u_z() * *x.get_u_z()
00615                   + f_u  * *x.get_u_zz()
00616                   + f_vv *  *y.get_u_z() * *y.get_u_z()
00617                   + f_v  * *y.get_u_zz()
00618                   + f_uv * (*x.get_u_z() * *y.get_u_z()
00619                          +  *x.get_u_z() * *y.get_u_z());
00620 
00621 
00622     *z.get_u_xxx() = f_uuu * cube(*x.get_u_x())
00623                    + 3.0 * f_uu * *x.get_u_x() * *x.get_u_xx()
00624                    + f_u * *x.get_u_xxx()
00625                    + f_vvv * cube(*y.get_u_x())
00626                    + 3.0 * f_vv * *y.get_u_x() * *y.get_u_xx()
00627                    + f_v * *y.get_u_xxx()
00628                    + 3.0 * f_uuv * square(*x.get_u_x()) * *y.get_u_x()
00629                    + 3.0 * f_uvv * *x.get_u_x()* square(*y.get_u_x())
00630                    + 3.0* f_uv * *x.get_u_xx() * *y.get_u_x()
00631                    + 3.0* f_uv * *x.get_u_x() * *y.get_u_xx();
00632 
00633     *z.get_u_xxy() = f_uuu * square(*x.get_u_x())* *x.get_u_y()
00634                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xy()
00635                    + f_uu * *x.get_u_y() * *x.get_u_xx()
00636                    + f_u * *x.get_u_xxy()
00637                    + f_vvv * square(*y.get_u_x())* *y.get_u_y()
00638                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xy()
00639                    + f_vv * *y.get_u_xx() * *y.get_u_y()
00640                    + f_v * *y.get_u_xxy()
00641                    + f_uuv * square(*x.get_u_x()) * *y.get_u_y()
00642                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_y() * *y.get_u_x()
00643                    + f_uvv * *x.get_u_y()* square(*y.get_u_x())
00644                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_y()
00645                    + f_uv * *x.get_u_xx() * *y.get_u_y()
00646                    + f_uv * *x.get_u_y() * *y.get_u_xx()
00647                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_x()
00648                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xy();
00649 
00650     *z.get_u_xxz() = f_uuu * square(*x.get_u_x())* *x.get_u_z()
00651                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xz()
00652                    + f_uu * *x.get_u_z() * *x.get_u_xx()
00653                    + f_u * *x.get_u_xxz()
00654                    + f_vvv * square(*y.get_u_x())* *y.get_u_z()
00655                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xz()
00656                    + f_vv * *y.get_u_xx() * *y.get_u_z()
00657                    + f_v * *y.get_u_xxz()
00658                    + f_uuv * square(*x.get_u_x()) * *y.get_u_z()
00659                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_z() * *y.get_u_x()
00660                    + f_uvv * *x.get_u_z()* square(*y.get_u_x())
00661                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_z()
00662                    + f_uv * *x.get_u_xx() * *y.get_u_z()
00663                    + f_uv * *x.get_u_z() * *y.get_u_xx()
00664                    + 2.0* f_uv * *x.get_u_xz() * *y.get_u_x()
00665                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xz();
00666 
00667     *z.get_u_xyy() = f_uuu * square(*x.get_u_y())* *x.get_u_x()
00668                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_xy()
00669                    + f_uu * *x.get_u_x() * *x.get_u_yy()
00670                    + f_u * *x.get_u_xyy()
00671                    + f_vvv * square(*y.get_u_y())* *y.get_u_x()
00672                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_xy()
00673                    + f_vv * *y.get_u_yy() * *y.get_u_x()
00674                    + f_v * *y.get_u_xyy()
00675                    + f_uuv * square(*x.get_u_y()) * *y.get_u_x()
00676                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_x() * *y.get_u_y()
00677                    + f_uvv * *x.get_u_x()* square(*y.get_u_y())
00678                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_x()
00679                    + f_uv * *x.get_u_yy() * *y.get_u_x()
00680                    + f_uv * *x.get_u_x() * *y.get_u_yy()
00681                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_y()
00682                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_xy();
00683 
00684     *z.get_u_xyz() = f_uuu * *x.get_u_y()* *x.get_u_x()* *x.get_u_z()
00685                    + f_uuv * *x.get_u_y() * *x.get_u_x()* *y.get_u_z()
00686                    + f_uu * *x.get_u_y() * *x.get_u_xz()
00687                    + f_uu * *x.get_u_x() * *x.get_u_yz()
00688 
00689                    + f_uuv * *y.get_u_y() * *x.get_u_x()* *x.get_u_z()
00690                    + f_uvv * *y.get_u_y() * *x.get_u_x()* *y.get_u_z()
00691                    + f_uv * *x.get_u_xz() * *y.get_u_y()
00692                    + f_uv * *y.get_u_yz() * *x.get_u_x()
00693 
00694                    + f_uu * *x.get_u_xy() * *x.get_u_z()
00695                    + f_uv * *x.get_u_xy() * *y.get_u_z()
00696                    + f_u * *x.get_u_xyz()
00697 
00698                    + f_uuv * *x.get_u_y() * *y.get_u_x()* *x.get_u_z()
00699                    + f_uvv * *x.get_u_y() * *y.get_u_x()* *y.get_u_z()
00700                    + f_uv * *y.get_u_x() * *x.get_u_yz()
00701                    + f_uv * *x.get_u_y() * *y.get_u_xz()
00702 
00703                    + f_uvv * *y.get_u_y() * *y.get_u_x()* *x.get_u_z()
00704                    + f_vvv * *y.get_u_y() * *y.get_u_x()* *y.get_u_z()
00705 
00706                    + f_vv * *y.get_u_xz() * *y.get_u_y()
00707                    + f_vv * *y.get_u_yz() * *y.get_u_x()
00708                    + f_vv * *y.get_u_xy() * *y.get_u_z()
00709                    + f_uv * *y.get_u_xy() * *x.get_u_z()
00710                    + f_v * *y.get_u_xyz();
00711 
00712     *z.get_u_xzz() = f_uuu * square(*x.get_u_z())* *x.get_u_x()
00713                    + 2.0 * f_uu * *x.get_u_z() * *x.get_u_xz()
00714                    + f_uu * *x.get_u_x() * *x.get_u_zz()
00715                    + f_u * *x.get_u_xzz()
00716                    + f_vvv * square(*y.get_u_z())* *y.get_u_x()
00717                    + 2.0 * f_vv * *y.get_u_z() * *y.get_u_xz()
00718                    + f_vv * *y.get_u_zz() * *y.get_u_x()
00719                    + f_v * *y.get_u_xzz()
00720                    + f_uuv * square(*x.get_u_z()) * *y.get_u_x()
00721                    + 2.0* f_uuv * *x.get_u_z() * *x.get_u_x() * *y.get_u_z()
00722                    + f_uvv * *x.get_u_x()* square(*y.get_u_z())
00723                    + 2.0 * f_uvv * *x.get_u_z()* *y.get_u_z() * *y.get_u_x()
00724                    + f_uv * *x.get_u_zz() * *y.get_u_x()
00725                    + f_uv * *x.get_u_x() * *y.get_u_zz()
00726                    + 2.0* f_uv * *x.get_u_xz() * *y.get_u_z()
00727                    + 2.0* f_uv * *x.get_u_z() * *y.get_u_xz();
00728 
00729     *z.get_u_yyy() = f_uuu * cube(*x.get_u_y())
00730                    + 3.0 * f_uu * *x.get_u_y() * *x.get_u_yy()
00731                    + f_u * *x.get_u_yyy()
00732                    + f_vvv * cube(*y.get_u_y())
00733                    + 3.0 * f_vv * *y.get_u_y() * *y.get_u_yy()
00734                    + f_v * *y.get_u_yyy()
00735                    + 3.0 * f_uuv * square(*x.get_u_y()) * *y.get_u_y()
00736                    + 3.0 * f_uvv * *x.get_u_y()* square(*y.get_u_y())
00737                    + 3.0 * f_uv * *x.get_u_yy() * *y.get_u_y()
00738                    + 3.0 * f_uv * *x.get_u_y() * *y.get_u_yy();
00739 
00740     *z.get_u_yyz() = f_uuu * square(*x.get_u_y())* *x.get_u_z()
00741                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_yz()
00742                    + f_uu * *x.get_u_z() * *x.get_u_yy()
00743                    + f_u * *x.get_u_yyz()
00744                    + f_vvv * square(*y.get_u_y())* *y.get_u_z()
00745                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_yz()
00746                    + f_vv * *y.get_u_yy() * *y.get_u_z()
00747                    + f_v * *y.get_u_yyz()
00748                    + f_uuv * square(*x.get_u_y()) * *y.get_u_z()
00749                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_z() * *y.get_u_y()
00750                    + f_uvv * *x.get_u_z()* square(*y.get_u_y())
00751                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_z()
00752                    + f_uv * *x.get_u_yy() * *y.get_u_z()
00753                    + f_uv * *x.get_u_z() * *y.get_u_yy()
00754                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_y()
00755                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_yz();
00756 
00757     *z.get_u_yyz() = f_uuu * square(*x.get_u_y())* *x.get_u_z()
00758                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_yz()
00759                    + f_uu * *x.get_u_z() * *x.get_u_yy()
00760                    + f_u * *x.get_u_yyz()
00761                    + f_vvv * square(*y.get_u_y())* *y.get_u_z()
00762                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_yz()
00763                    + f_vv * *y.get_u_yy() * *y.get_u_z()
00764                    + f_v * *y.get_u_yyz()
00765                    + f_uuv * square(*x.get_u_y()) * *y.get_u_z()
00766                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_z() * *y.get_u_y()
00767                    + f_uvv * *x.get_u_z()* square(*y.get_u_y())
00768                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_z()
00769                    + f_uv * *x.get_u_yy() * *y.get_u_z()
00770                    + f_uv * *x.get_u_z() * *y.get_u_yy()
00771                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_y()
00772                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_yz();
00773 
00774     *z.get_u_yzz() = f_uuu * square(*x.get_u_z())* *x.get_u_y()
00775                    + 2.0 * f_uu * *x.get_u_z() * *x.get_u_yz()
00776                    + f_uu * *x.get_u_y() * *x.get_u_zz()
00777                    + f_u * *x.get_u_yzz()
00778                    + f_vvv * square(*y.get_u_z())* *y.get_u_y()
00779                    + 2.0 * f_vv * *y.get_u_z() * *y.get_u_yz()
00780                    + f_vv * *y.get_u_zz() * *y.get_u_y()
00781                    + f_v * *y.get_u_yzz()
00782                    + f_uuv * square(*x.get_u_z()) * *y.get_u_y()
00783                    + 2.0* f_uuv * *x.get_u_z() * *x.get_u_y() * *y.get_u_z()
00784                    + f_uvv * *x.get_u_y()* square(*y.get_u_z())
00785                    + 2.0 * f_uvv * *x.get_u_z()* *y.get_u_z() * *y.get_u_y()
00786                    + f_uv * *x.get_u_zz() * *y.get_u_y()
00787                    + f_uv * *x.get_u_y() * *y.get_u_zz()
00788                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_z()
00789                    + 2.0* f_uv * *x.get_u_z() * *y.get_u_yz();
00790 
00791     *z.get_u_zzz() = f_uuu * cube(*x.get_u_z())
00792                    + 3.0 * f_uu * *x.get_u_z() * *x.get_u_zz()
00793                    + f_u * *x.get_u_zzz()
00794                    + f_vvv * cube(*y.get_u_z())
00795                    + 3.0 * f_vv * *y.get_u_z() * *y.get_u_zz()
00796                    + f_v * *y.get_u_zzz()
00797                    + 3.0 * f_uuv * square(*x.get_u_z()) * *y.get_u_z()
00798                    + 3.0 * f_uvv * *x.get_u_z()* square(*y.get_u_z())
00799                    + 3.0 * f_uv * *x.get_u_zz() * *y.get_u_z()
00800                    + 3.0 * f_uv * *x.get_u_z() * *y.get_u_zz();
00801 }
00802 
00807   df3_three_variable sqrt(const df3_three_variable& x)
00808   {
00809     df3_three_variable z;
00810     double u=::sqrt(*x.get_u());
00811     *z.get_u()=u;
00812     double xinv=1.0/(*x.get_u());
00813     double zp=0.5/u;
00814     double zp2=-0.5*zp*xinv;
00815     double zp3=-1.5*zp2*xinv;
00816 
00817 
00818     set_derivatives(z,x,u,zp,zp2,zp3);
00819 
00820     return z;
00821   }
00822 
00827   df3_three_variable atan(const df3_three_variable& x)
00828   {
00829     df3_three_variable z;
00830     double cx=value(x);
00831     double d=1.0/(1+square(cx));
00832     double d2=square(d);
00833     double u=::atan(cx);
00834     *z.get_u()=u;
00835     double zp=d;
00836     double zp2=-2.0*cx*d2;
00837     double zp3=-2.0*d2+8*cx*cx*d*d2;
00838 
00839     set_derivatives(z,x,u,zp,zp2,zp3);
00840     return z;
00841   }
00842 
00847   df3_three_variable pow(const df3_three_variable& x,
00848     const df3_three_variable& y)
00849   {
00850     df3_three_variable z;
00851     z=exp(y*log(x));
00852     return z;
00853   }
00854 
00859   df3_three_variable square(const df3_three_variable& x)
00860   {
00861     df3_three_variable z;
00862     double u=value(x);
00863     *z.get_u()=u*u;
00864     double zp=2.0*u;
00865     double zp2=2.0;
00866     double zp3=0.0;
00867 
00868     set_derivatives(z,x,u,zp,zp2,zp3);
00869     return z;
00870   }
00871 
00876   df3_three_variable tan(const df3_three_variable& x)
00877   {
00878     df3_three_variable z;
00879     double u=::tan(*x.get_u());
00880     *z.get_u()=u;
00881     double v=1.0/::cos(*x.get_u());
00882     double w=::sin(*x.get_u());
00883     double v2=v*v;
00884     double zp=v2;
00885     double zp2=2.0*w*v2*v;
00886     double zp3=(4.0*w*w+2.0)*v2*v2;
00887 
00888     set_derivatives(z,x,u,zp,zp2,zp3);
00889     return z;
00890   }
00891 
00896   df3_three_variable sin(const df3_three_variable& x)
00897   {
00898     df3_three_variable z;
00899     double u=::sin(*x.get_u());
00900     *z.get_u()=u;
00901     double zp=::cos(*x.get_u());
00902     double zp2=-u;
00903     double zp3=-zp;
00904 
00905     set_derivatives(z,x,u,zp,zp2,zp3);
00906     return z;
00907   }
00908 
00913   df3_three_variable log(const df3_three_variable& x)
00914   {
00915     df3_three_variable z;
00916     double u=::log(*x.get_u());
00917     *z.get_u()=u;
00918     double zp=1/(*x.get_u());
00919     double zp2=-zp*zp;
00920     double zp3=-2.0*zp*zp2;
00921 
00922     set_derivatives(z,x,u,zp,zp2,zp3);
00923     return z;
00924   }
00925 
00930   df3_three_variable exp(const df3_three_variable& x)
00931   {
00932     df3_three_variable z;
00933     double u=::exp(*x.get_u());
00934     *z.get_u()=u;
00935     double zp=u;
00936     double zp2=u;
00937     double zp3=u;
00938 
00939     set_derivatives(z,x,u,zp,zp2,zp3);
00940     return z;
00941   }
00942 
00947   df3_three_variable inv(const df3_three_variable& x)
00948   {
00949     df3_three_variable z;
00950     double xinv=1.0/(*x.get_u());
00951     *z.get_u()=xinv;
00952     double zp=-xinv*xinv;
00953     double zp2=-2.0*zp*xinv;
00954     double zp3=-3.0*zp2*xinv;
00955 
00956     set_derivatives(z,x,xinv,zp,zp2,zp3);
00957 
00958     return z;
00959   }
00963 df3_three_variable& df3_three_variable::operator=(const df3_three_variable& _v)
00964 {
00965   *get_u() =  *_v.get_u();
00966   *get_u_x() = *_v.get_u_x();
00967   *get_u_y() = *_v.get_u_y();
00968   *get_u_z() = *_v.get_u_z();
00969   *get_u_xx() = *_v.get_u_xx();
00970   *get_u_xy() = *_v.get_u_xy();
00971   *get_u_xz() = *_v.get_u_xz();
00972   *get_u_yy() = *_v.get_u_yy();
00973   *get_u_yz() = *_v.get_u_yz();
00974   *get_u_zz() = *_v.get_u_zz();
00975   *get_u_xxx() = *_v.get_u_xxx();
00976   *get_u_xxy() = *_v.get_u_xxy();
00977   *get_u_xxz() = *_v.get_u_xxz();
00978   *get_u_xyy() = *_v.get_u_xyy();
00979   *get_u_xyz() = *_v.get_u_xyz();
00980   *get_u_xzz() = *_v.get_u_xzz();
00981   *get_u_yyy() = *_v.get_u_yyy();
00982   *get_u_yyz() = *_v.get_u_yyz();
00983   *get_u_yzz() = *_v.get_u_yzz();
00984   *get_u_zzz() = *_v.get_u_zzz();
00985 
00986   return *this;
00987 }
00991 df3_three_variable& df3_three_variable::operator=(double x)
00992 {
00993   *get_u() = x;
00994   *get_u_x() = 0.0;
00995   *get_u_y() = 0.0;
00996   *get_u_z() = 0.0;
00997   *get_u_xx() = 0.0;
00998   *get_u_xy() = 0.0;
00999   *get_u_xz() = 0.0;
01000   *get_u_yy() = 0.0;
01001   *get_u_yz() = 0.0;
01002   *get_u_zz() = 0.0;
01003   *get_u_xxx() = 0.0;
01004   *get_u_xxy() = 0.0;
01005   *get_u_xxz() = 0.0;
01006   *get_u_xyy() = 0.0;
01007   *get_u_xyz() = 0.0;
01008   *get_u_xzz() = 0.0;
01009   *get_u_yyy() = 0.0;
01010   *get_u_yyz() = 0.0;
01011   *get_u_yzz() = 0.0;
01012   *get_u_zzz() = 0.0;
01013 
01014   return *this;
01015 }
01019 df3_three_variable operator*(
01020   const df3_three_variable& x,
01021   const df3_three_variable& y)
01022 {
01023   df3_three_variable z;
01024   double u = *x.get_u() * *y.get_u();
01025   *z.get_u() = u;
01026   double f_u = *y.get_u();
01027   double f_v = *x.get_u();
01028   double f_uu = 0.0;
01029   double f_uv = 1.0;
01030   double f_vv = 0.0;
01031   double f_uuu = 0.0;
01032   double f_uuv = 0.0;
01033   double f_uvv = 0.0;
01034   double f_vvv = 0.0;
01035   set_derivatives(z,x,y,u,
01036     f_u, f_v,
01037     f_uu, f_uv, f_vv,
01038     f_uuu, f_uuv, f_uvv, f_vvv);
01039   return z;
01040 }
01041 
01046   df3_three_variable operator * (double x,
01047     const df3_three_variable& v)
01048   {
01049     df3_three_variable z;
01050 
01051     *z.get_u() = x * *v.get_u();
01052     *z.get_u_x() = x * *v.get_u_x();
01053     *z.get_u_y() = x * *v.get_u_y();
01054     *z.get_u_z() = x * *v.get_u_z();
01055     *z.get_u_xx() = x * *v.get_u_xx();
01056     *z.get_u_xy() = x * *v.get_u_xy();
01057     *z.get_u_xz() = x * *v.get_u_xz();
01058     *z.get_u_yy() = x * *v.get_u_yy();
01059     *z.get_u_yz() = x * *v.get_u_yz();
01060     *z.get_u_zz() = x * *v.get_u_zz();
01061     *z.get_u_xxx() = x * *v.get_u_xxx();
01062     *z.get_u_xxy() = x * *v.get_u_xxy();
01063     *z.get_u_xxz() = x * *v.get_u_xxz();
01064     *z.get_u_xyy() = x * *v.get_u_xyy();
01065     *z.get_u_xyz() = x * *v.get_u_xyz();
01066     *z.get_u_xzz() = x * *v.get_u_xzz();
01067     *z.get_u_yyy() = x * *v.get_u_yyy();
01068     *z.get_u_yyz() = x * *v.get_u_yyz();
01069     *z.get_u_yzz() = x * *v.get_u_yzz();
01070     *z.get_u_zzz() = x * *v.get_u_zzz();
01071 
01072     return z;
01073   }
01074 
01079   df3_three_variable fabs(const df3_three_variable& v)
01080   {
01081     df3_three_variable z;
01082     if (value(v)>=0)
01083     {
01084       *z.get_u() = *v.get_u();
01085       *z.get_u_x() = *v.get_u_x();
01086       *z.get_u_y() = *v.get_u_y();
01087       *z.get_u_z() = *v.get_u_z();
01088       *z.get_u_xx() = *v.get_u_xx();
01089       *z.get_u_xy() = *v.get_u_xy();
01090       *z.get_u_xz() = *v.get_u_xz();
01091       *z.get_u_yy() = *v.get_u_yy();
01092       *z.get_u_yz() = *v.get_u_yz();
01093       *z.get_u_zz() = *v.get_u_zz();
01094       *z.get_u_xxx() = *v.get_u_xxx();
01095       *z.get_u_xxy() = *v.get_u_xxy();
01096       *z.get_u_xxz() = *v.get_u_xxz();
01097       *z.get_u_xyy() = *v.get_u_xyy();
01098       *z.get_u_xyz() = *v.get_u_xyz();
01099       *z.get_u_xzz() = *v.get_u_xzz();
01100       *z.get_u_yyy() = *v.get_u_yyy();
01101       *z.get_u_yyz() = *v.get_u_yyz();
01102       *z.get_u_yzz() = *v.get_u_yzz();
01103       *z.get_u_zzz() = *v.get_u_zzz();
01104     }
01105     else
01106     {
01107       *z.get_u() = -*v.get_u();
01108       *z.get_u_x() = -*v.get_u_x();
01109       *z.get_u_y() = -*v.get_u_y();
01110       *z.get_u_z() = -*v.get_u_z();
01111       *z.get_u_xx() = -*v.get_u_xx();
01112       *z.get_u_xy() = -*v.get_u_xy();
01113       *z.get_u_xz() = -*v.get_u_xz();
01114       *z.get_u_yy() = -*v.get_u_yy();
01115       *z.get_u_yz() = -*v.get_u_yz();
01116       *z.get_u_zz() = -*v.get_u_zz();
01117       *z.get_u_xxx() = -*v.get_u_xxx();
01118       *z.get_u_xxy() = -*v.get_u_xxy();
01119       *z.get_u_xxz() = -*v.get_u_xxz();
01120       *z.get_u_xyy() = -*v.get_u_xyy();
01121       *z.get_u_xyz() = -*v.get_u_xyz();
01122       *z.get_u_xzz() = -*v.get_u_xzz();
01123       *z.get_u_yyy() = -*v.get_u_yyy();
01124       *z.get_u_yyz() = -*v.get_u_yyz();
01125       *z.get_u_yzz() = -*v.get_u_yzz();
01126       *z.get_u_zzz() = -*v.get_u_zzz();
01127     }
01128 
01129     return z;
01130   }
01131 
01136   df3_three_variable operator * (const df3_three_variable& v,
01137     double x)
01138   {
01139     df3_three_variable z;
01140 
01141     *z.get_u() = x * *v.get_u();
01142     *z.get_u_x() = x * *v.get_u_x();
01143     *z.get_u_y() = x * *v.get_u_y();
01144     *z.get_u_z() = x * *v.get_u_z();
01145     *z.get_u_xx() = x * *v.get_u_xx();
01146     *z.get_u_xy() = x * *v.get_u_xy();
01147     *z.get_u_xz() = x * *v.get_u_xz();
01148     *z.get_u_yy() = x * *v.get_u_yy();
01149     *z.get_u_yz() = x * *v.get_u_yz();
01150     *z.get_u_zz() = x * *v.get_u_zz();
01151     *z.get_u_xxx() = x * *v.get_u_xxx();
01152     *z.get_u_xxy() = x * *v.get_u_xxy();
01153     *z.get_u_xxz() = x * *v.get_u_xxz();
01154     *z.get_u_xyy() = x * *v.get_u_xyy();
01155     *z.get_u_xyz() = x * *v.get_u_xyz();
01156     *z.get_u_xzz() = x * *v.get_u_xzz();
01157     *z.get_u_yyy() = x * *v.get_u_yyy();
01158     *z.get_u_yyz() = x * *v.get_u_yyz();
01159     *z.get_u_yzz() = x * *v.get_u_yzz();
01160     *z.get_u_zzz() = x * *v.get_u_zzz();
01161 
01162     return z;
01163   }
01164 
01169   df3_three_variable operator / (const df3_three_variable& x,
01170     double y)
01171   {
01172     double u=1/y;
01173     return x*u;
01174   }
01175 
01180   df3_three_variable operator / (const df3_three_variable& x,
01181     const df3_three_variable& y)
01182   {
01183     df3_three_variable u=inv(y);
01184     return x*u;
01185   }
01186 
01191   df3_three_variable operator / (const double x,
01192     const df3_three_variable& y)
01193   {
01194     df3_three_variable u=inv(y);
01195     return x*u;
01196   }
01197 
01198  /*
01199   df3_three_variable operator / (const df3_three_variable& x,
01200     const df3_three_variable& y)
01201   {
01202     df3_three_variable z;
01203     double yinv =  1.0 / (*y.get_u());
01204     double yinv2 = yinv * yinv;
01205     double yinv3 = yinv * yinv2;
01206     doubl yd = *y.get_udot();
01207 
01208     // *z.get_u() = *x.get_u() /  *y.get_u();
01209     *z.get_u() = *x.get_u() * yinv;
01210 
01211     *z.get_udot() =  - (*x.get_u()) * yinv2 * yd
01212                   + *x.get_udot() * yinv;
01213 
01214     *z.get_udot2() = *x.get_udot2() * yinv
01215                    - 2.0 * *x.get_udot() * yd * yinv2
01216                    + 2.0 * *x.get_u() * yinv3  * yd *yd
01217                    -  *x.get_u() * yinv2 * y.get_udot2();
01218 
01219     *z.get_udot3() = *x.get_udot3() * yinv
01220                    + 3.0 * *x.get_udot2() * *y.get_udot()
01221                    + 3.0 * *x.get_udot() * *y.get_udot2()
01222                    +  *x.get_u() * *y.get_udot3();
01223   }
01224  */
01225 
01230   df3_three_variable operator + (const double x,const df3_three_variable& v)
01231   {
01232     df3_three_variable z;
01233     *z.get_u() = x + *v.get_u();
01234     *z.get_u_x() = *v.get_u_x();
01235     *z.get_u_y() = *v.get_u_y();
01236     *z.get_u_z() = *v.get_u_z();
01237     *z.get_u_xx() = *v.get_u_xx();
01238     *z.get_u_xy() = *v.get_u_xy();
01239     *z.get_u_xz() = *v.get_u_xz();
01240     *z.get_u_yy() = *v.get_u_yy();
01241     *z.get_u_yz() = *v.get_u_yz();
01242     *z.get_u_zz() = *v.get_u_zz();
01243     *z.get_u_xxx() = *v.get_u_xxx();
01244     *z.get_u_xxy() = *v.get_u_xxy();
01245     *z.get_u_xxz() = *v.get_u_xxz();
01246     *z.get_u_xyy() = *v.get_u_xyy();
01247     *z.get_u_xyz() = *v.get_u_xyz();
01248     *z.get_u_xzz() = *v.get_u_xzz();
01249     *z.get_u_yyy() = *v.get_u_yyy();
01250     *z.get_u_yyz() = *v.get_u_yyz();
01251     *z.get_u_yzz() = *v.get_u_yzz();
01252     *z.get_u_zzz() = *v.get_u_zzz();
01253 
01254     return z;
01255   }
01256 
01261   df3_three_variable operator + (const df3_three_variable& v,const double x)
01262   {
01263     df3_three_variable z;
01264 
01265     *z.get_u() = x + *v.get_u();
01266     *z.get_u_x() = *v.get_u_x();
01267     *z.get_u_y() = *v.get_u_y();
01268     *z.get_u_z() = *v.get_u_z();
01269     *z.get_u_xx() = *v.get_u_xx();
01270     *z.get_u_xy() = *v.get_u_xy();
01271     *z.get_u_xz() = *v.get_u_xz();
01272     *z.get_u_yy() = *v.get_u_yy();
01273     *z.get_u_yz() = *v.get_u_yz();
01274     *z.get_u_zz() = *v.get_u_zz();
01275     *z.get_u_xxx() = *v.get_u_xxx();
01276     *z.get_u_xxy() = *v.get_u_xxy();
01277     *z.get_u_xxz() = *v.get_u_xxz();
01278     *z.get_u_xyy() = *v.get_u_xyy();
01279     *z.get_u_xyz() = *v.get_u_xyz();
01280     *z.get_u_xzz() = *v.get_u_xzz();
01281     *z.get_u_yyy() = *v.get_u_yyy();
01282     *z.get_u_yyz() = *v.get_u_yyz();
01283     *z.get_u_yzz() = *v.get_u_yzz();
01284     *z.get_u_zzz() = *v.get_u_zzz();
01285 
01286 
01287 
01288     return z;
01289   }
01290 
01295   df3_three_variable operator + (const df3_three_variable& x,
01296     const df3_three_variable& y)
01297   {
01298     df3_three_variable z;
01299     *z.get_u() = *x.get_u() + *y.get_u();
01300     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
01301     *z.get_u_z() = *x.get_u_z() + *y.get_u_z();
01302     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
01303     *z.get_u_xx() = *x.get_u_xx()+ *y.get_u_xx();
01304     *z.get_u_xy() = *x.get_u_xy()+ *y.get_u_xy();
01305     *z.get_u_xz() = *x.get_u_xz()+ *y.get_u_xz();
01306     *z.get_u_yy() = *x.get_u_yy()+ *y.get_u_yy();
01307     *z.get_u_yz() = *x.get_u_yz()+ *y.get_u_yz();
01308     *z.get_u_zz() = *x.get_u_zz()+ *y.get_u_zz();
01309     *z.get_u_xxx() = *x.get_u_xxx()+ *y.get_u_xxx();
01310     *z.get_u_xxy() = *x.get_u_xxy()+ *y.get_u_xxy();
01311     *z.get_u_xxz() = *x.get_u_xxz()+ *y.get_u_xxz();
01312     *z.get_u_xyy() = *x.get_u_xyy()+ *y.get_u_xyy();
01313     *z.get_u_xyz() = *x.get_u_xyz()+ *y.get_u_xyz();
01314     *z.get_u_xzz() = *x.get_u_xzz()+ *y.get_u_xzz();
01315     *z.get_u_yyy() = *x.get_u_yyy()+ *y.get_u_yyy();
01316     *z.get_u_yyz() = *x.get_u_yyz()+ *y.get_u_yyz();
01317     *z.get_u_yzz() = *x.get_u_yzz()+ *y.get_u_yzz();
01318     *z.get_u_zzz() = *x.get_u_zzz()+ *y.get_u_zzz();
01319     return z;
01320   }
01321 
01326   df3_three_variable operator - (const df3_three_variable& x,
01327     const df3_three_variable& y)
01328   {
01329     df3_three_variable z;
01330     *z.get_u() = *x.get_u() - *y.get_u();
01331     *z.get_u_x() = *x.get_u_x() - *y.get_u_x();
01332     *z.get_u_z() = *x.get_u_z() - *y.get_u_z();
01333     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
01334     *z.get_u_xx() = *x.get_u_xx() - *y.get_u_xx();
01335     *z.get_u_xy() = *x.get_u_xy() - *y.get_u_xy();
01336     *z.get_u_xz() = *x.get_u_xz() - *y.get_u_xz();
01337     *z.get_u_yy() = *x.get_u_yy() - *y.get_u_yy();
01338     *z.get_u_yz() = *x.get_u_yz() - *y.get_u_yz();
01339     *z.get_u_zz() = *x.get_u_zz() - *y.get_u_zz();
01340     *z.get_u_xxx() = *x.get_u_xxx() - *y.get_u_xxx();
01341     *z.get_u_xxy() = *x.get_u_xxy() - *y.get_u_xxy();
01342     *z.get_u_xxz() = *x.get_u_xxz() - *y.get_u_xxz();
01343     *z.get_u_xyy() = *x.get_u_xyy() - *y.get_u_xyy();
01344     *z.get_u_xyz() = *x.get_u_xyz() - *y.get_u_xyz();
01345     *z.get_u_xzz() = *x.get_u_xzz() - *y.get_u_xzz();
01346     *z.get_u_yyy() = *x.get_u_yyy() - *y.get_u_yyy();
01347     *z.get_u_yyz() = *x.get_u_yyz() - *y.get_u_yyz();
01348     *z.get_u_yzz() = *x.get_u_yzz() - *y.get_u_yzz();
01349     *z.get_u_zzz() = *x.get_u_zzz() - *y.get_u_zzz();
01350     return z;
01351   }
01352 
01357   df3_three_variable operator - (double x,
01358     const df3_three_variable& v)
01359   {
01360     df3_three_variable z;
01361 
01362     *z.get_u() = x - *v.get_u();
01363     *z.get_u_x() = -*v.get_u_x();
01364     *z.get_u_y() = -*v.get_u_y();
01365     *z.get_u_z() = -*v.get_u_z();
01366     *z.get_u_xx() = -*v.get_u_xx();
01367     *z.get_u_xy() = -*v.get_u_xy();
01368     *z.get_u_xz() = -*v.get_u_xz();
01369     *z.get_u_yy() = -*v.get_u_yy();
01370     *z.get_u_yz() = -*v.get_u_yz();
01371     *z.get_u_zz() = -*v.get_u_zz();
01372     *z.get_u_xxx() = -*v.get_u_xxx();
01373     *z.get_u_xxy() = -*v.get_u_xxy();
01374     *z.get_u_xxz() = -*v.get_u_xxz();
01375     *z.get_u_xyy() = -*v.get_u_xyy();
01376     *z.get_u_xyz() = -*v.get_u_xyz();
01377     *z.get_u_xzz() = -*v.get_u_xzz();
01378     *z.get_u_yyy() = -*v.get_u_yyy();
01379     *z.get_u_yyz() = -*v.get_u_yyz();
01380     *z.get_u_yzz() = -*v.get_u_yzz();
01381     *z.get_u_zzz() = -*v.get_u_zzz();
01382 
01383     return z;
01384   }
01385 
01390   df3_three_variable operator - (const df3_three_variable& v,
01391     double y)
01392   {
01393     df3_three_variable z;
01394 
01395     *z.get_u() =  *v.get_u()-y;
01396     *z.get_u_x() = *v.get_u_x();
01397     *z.get_u_y() = *v.get_u_y();
01398     *z.get_u_z() = *v.get_u_z();
01399     *z.get_u_xx() = *v.get_u_xx();
01400     *z.get_u_xy() = *v.get_u_xy();
01401     *z.get_u_xz() = *v.get_u_xz();
01402     *z.get_u_yy() = *v.get_u_yy();
01403     *z.get_u_yz() = *v.get_u_yz();
01404     *z.get_u_zz() = *v.get_u_zz();
01405     *z.get_u_xxx() = *v.get_u_xxx();
01406     *z.get_u_xxy() = *v.get_u_xxy();
01407     *z.get_u_xxz() = *v.get_u_xxz();
01408     *z.get_u_xyy() = *v.get_u_xyy();
01409     *z.get_u_xyz() = *v.get_u_xyz();
01410     *z.get_u_xzz() = *v.get_u_xzz();
01411     *z.get_u_yyy() = *v.get_u_yyy();
01412     *z.get_u_yyz() = *v.get_u_yyz();
01413     *z.get_u_yzz() = *v.get_u_yzz();
01414     *z.get_u_zzz() = *v.get_u_zzz();
01415 
01416     return z;
01417   }
01421 init_df3_three_variable::~init_df3_three_variable()
01422 {
01423   num_local_ind_var--;
01424   if (num_local_ind_var<0)
01425   {
01426     cerr << "This can't happen" << endl;
01427     ad_exit(1);
01428   }
01429   if (num_local_ind_var==0)
01430   {
01431     num_ind_var=0;
01432   }
01433 }
01437   init_df3_three_variable::init_df3_three_variable(const df1b2variable& _v)
01438   {
01439     ADUNCONST(df1b2variable,v)
01440     if (num_local_ind_var>2)
01441     {
01442       cerr << "can only have 3 independent_variables in df3_three_variable"
01443        " function" << endl;
01444       ad_exit(1);
01445     }
01446     ind_var[num_ind_var++]=&v;
01447     num_local_ind_var++;
01448     *get_u() =  *v.get_u();
01449     *get_u_x() = 0.0;
01450     *get_u_y() = 0.0;
01451     *get_u_z() = 0.0;
01452     switch(num_local_ind_var)
01453     {
01454     case 1:
01455       *get_u_x() = 1.0;
01456       break;
01457     case 2:
01458       *get_u_y() = 1.0;
01459       break;
01460     case 3:
01461       *get_u_z() = 1.0;
01462       break;
01463     default:
01464       cerr << "illegal num_ind_var value of " << num_ind_var
01465            << " in  df3_three_variable function" << endl;
01466       ad_exit(1);
01467     }
01468     *get_u_xx() = 0.0;
01469     *get_u_xy() = 0.0;
01470     *get_u_xz() = 0.0;
01471     *get_u_yy() = 0.0;
01472     *get_u_yz() = 0.0;
01473     *get_u_zz() = 0.0;
01474     *get_u_xxx() = 0.0;
01475     *get_u_xxy() = 0.0;
01476     *get_u_xxz() = 0.0;
01477     *get_u_xyy() = 0.0;
01478     *get_u_xyz() = 0.0;
01479     *get_u_xzz() = 0.0;
01480     *get_u_yyy() = 0.0;
01481     *get_u_yyz() = 0.0;
01482     *get_u_yzz() = 0.0;
01483     *get_u_zzz() = 0.0;
01484   }
01485 
01490   init_df3_three_variable::init_df3_three_variable(double v)
01491   {
01492     if (num_local_ind_var>2)
01493     {
01494       cerr << "can only have 3 independent_variables in df3_three_variable"
01495        " function" << endl;
01496       ad_exit(1);
01497     }
01498     num_local_ind_var++;
01499     *get_u() =  v;
01500     *get_u_x() = 0.0;
01501     *get_u_y() = 0.0;
01502     *get_u_z() = 0.0;
01503     switch(num_local_ind_var)
01504     {
01505     case 1:
01506       *get_u_x() = 1.0;
01507       break;
01508     case 2:
01509       *get_u_y() = 1.0;
01510       break;
01511     case 3:
01512       *get_u_z() = 1.0;
01513       break;
01514     default:
01515       cerr << "illegal num_ind_var value of " << num_ind_var
01516            << " in  df3_three_variable function" << endl;
01517       ad_exit(1);
01518     }
01519     *get_u_xx() = 0.0;
01520     *get_u_xy() = 0.0;
01521     *get_u_xz() = 0.0;
01522     *get_u_yy() = 0.0;
01523     *get_u_yz() = 0.0;
01524     *get_u_zz() = 0.0;
01525     *get_u_xxx() = 0.0;
01526     *get_u_xxy() = 0.0;
01527     *get_u_xxz() = 0.0;
01528     *get_u_xyy() = 0.0;
01529     *get_u_xyz() = 0.0;
01530     *get_u_xzz() = 0.0;
01531     *get_u_yyy() = 0.0;
01532     *get_u_yyz() = 0.0;
01533     *get_u_yzz() = 0.0;
01534     *get_u_zzz() = 0.0;
01535   }
01540 df3_three_matrix choleski_decomp(const df3_three_matrix& MM)
01541 {
01542   // kludge to deal with constantness
01543   df3_three_matrix & M= (df3_three_matrix &) MM;
01544   int rmin=M.indexmin();
01545   int cmin=M(rmin).indexmin();
01546   int rmax=M.indexmax();
01547   int cmax=M(rmin).indexmax();
01548   if (rmin !=1 || cmin !=1)
01549   {
01550     cerr << "minimum row and column inidices must equal 1 in "
01551       "df1b2matrix choleski_decomp(const df3_three_atrix& MM)"
01552          << endl;
01553     ad_exit(1);
01554   }
01555   if (rmax !=cmax)
01556   {
01557     cerr << "Error in df1b2matrix choleski_decomp(const df3_three_matrix& MM)"
01558       " Matrix not square" << endl;
01559     ad_exit(1);
01560   }
01561 
01562   int n=rmax-rmin+1;
01563   df3_three_matrix L(1,n,1,n);
01564 #ifndef SAFE_INITIALIZE
01565     L.initialize();
01566 #endif
01567 
01568   int i,j,k;
01569   df3_three_variable tmp;
01570 
01571     if (value(M(1,1))<=0)
01572     {
01573       cerr << "Error matrix not positive definite in choleski_decomp"
01574         <<endl;
01575       ad_exit(1);
01576     }
01577 
01578   L(1,1)=sqrt(M(1,1));
01579   for (i=2;i<=n;i++)
01580   {
01581     L(i,1)=M(i,1)/L(1,1);
01582   }
01583 
01584   for (i=2;i<=n;i++)
01585   {
01586     for (j=2;j<=i-1;j++)
01587     {
01588       tmp=M(i,j);
01589       for (k=1;k<=j-1;k++)
01590       {
01591         tmp-=L(i,k)*L(j,k);
01592       }
01593       L(i,j)=tmp/L(j,j);
01594     }
01595     tmp=M(i,i);
01596     for (k=1;k<=i-1;k++)
01597     {
01598       tmp-=L(i,k)*L(i,k);
01599     }
01600 
01601     if (value(tmp)<=0)
01602     {
01603       cerr << "Error matrix not positive definite in choleski_decomp"
01604         <<endl;
01605       ad_exit(1);
01606     }
01607 
01608     L(i,i)=sqrt(tmp);
01609   }
01610 
01611   return L;
01612 }