ADMB Documentation  11.1.1914
 All Classes Files Functions Variables Typedefs Friends Defines
df33fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df33fun.cpp 1709 2014-02-28 21:48:21Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00011 #include <df1b2fun.h>
00012 #include "df33fun.h"
00013   df1b2variable * df3_three_variable::ind_var[3];
00014   int df3_three_variable::num_ind_var=0;
00015   int df3_three_variable::num_local_ind_var=0;
00016 
00021   df3_three_variable::df3_three_variable(const df3_three_variable& x)
00022   {
00023     memcpy(&(v[0]),&(x.v[0]),20*sizeof(double));
00024   }
00025 
00030  df3_three_vector::df3_three_vector(const df3_three_vector& m2)
00031  {
00032    index_min=m2.index_min;
00033    index_max=m2.index_max;
00034    shape=m2.shape;
00035    if (shape)
00036    {
00037      (shape->ncopies)++;
00038    }
00039    v = m2.v;
00040  }
00041 
00046 df1b2variable& df1b2variable::operator = (const df3_three_variable& v)
00047 {
00048   df1b2variable * px=df3_three_variable::ind_var[0];
00049   df1b2variable * py=df3_three_variable::ind_var[1];
00050   df1b2variable * pv=df3_three_variable::ind_var[2];
00051   df3_three_variable::num_ind_var=0;
00052   df3_three_variable::ind_var[0]=0;
00053   df3_three_variable::ind_var[1]=0;
00054   df3_three_variable::ind_var[2]=0;
00055   //df1b2variable * px=0;
00056   double  dfx= *v.get_u_x();
00057   double  dfy= *v.get_u_y();
00058   double  dfz= *v.get_u_z();
00059   double * xd=px->get_u_dot();
00060   double * yd=py->get_u_dot();
00061   double * vd=pv->get_u_dot();
00062   double * zd=get_u_dot();
00063   *get_u()=*v.get_u();
00064   for (int i=0;i<df1b2variable::nvar;i++)
00065   {
00066     *zd++ = dfx * *xd++ + dfy * *yd++ + dfz * *vd++;
00067   }
00068 
00069  /*
00070   cout << *v.get_u()  << " ";
00071   cout << *v.get_udot()  << " ";
00072   cout << *v.get_udot2()  << " ";
00073   cout << *v.get_udot3()  << endl;
00074   */
00075   f1b2gradlist->write_pass1(px,
00076     py,
00077     pv,
00078     this,
00079     *v.get_u_x(),
00080     *v.get_u_y(),
00081     *v.get_u_z(),
00082     *v.get_u_xx(),
00083     *v.get_u_xy(),
00084     *v.get_u_xz(),
00085     *v.get_u_yy(),
00086     *v.get_u_yz(),
00087     *v.get_u_zz(),
00088     *v.get_u_xxx(),
00089     *v.get_u_xxy(),
00090     *v.get_u_xxz(),
00091     *v.get_u_xyy(),
00092     *v.get_u_xyz(),
00093     *v.get_u_xzz(),
00094     *v.get_u_yyy(),
00095     *v.get_u_yyz(),
00096     *v.get_u_yzz(),
00097     *v.get_u_zzz());
00098   return *this;
00099 }
00100 
00105 void df3_three_variable::initialize(void)
00106 {
00107   for (int i=0;i<20;i++)
00108     v[i]=0.0;
00109 }
00110 
00115  df3_three_vector::~df3_three_vector()
00116  {
00117    deallocate();
00118  }
00119 
00124  void df3_three_vector::deallocate(void)
00125  {
00126    if(shape)
00127    {
00128      if (shape->ncopies)
00129      {
00130        (shape->ncopies)--;
00131      }
00132      else
00133      {
00134        v = (df3_three_variable*) (shape->trueptr);
00135        delete [] v;
00136        v = NULL;
00137        delete shape;
00138        shape=0;
00139      }
00140    }
00141  }
00142 
00147  dvector value(const df3_three_vector& v)
00148  {
00149    int mmin=v.indexmin();
00150    int mmax=v.indexmax();
00151    dvector cv(mmin,mmax);
00152    for (int i=mmin;i<=mmax;i++)
00153    {
00154      cv(i)=value(v(i));
00155    }
00156    return cv;
00157  }
00158 
00163   void df3_three_vector::initialize(void)
00164   {
00165     int mmin=indexmin();
00166     int mmax=indexmax();
00167     for (int i=mmin;i<=mmax;i++)
00168     {
00169       (*this)(i)=0.0;
00170     }
00171   }
00172 
00177   df3_three_vector::df3_three_vector(void)
00178   {
00179     allocate();
00180   }
00181 
00186   df3_three_vector::df3_three_vector(int min,int max)
00187   {
00188     allocate(min,max);
00189   }
00190 
00195   void df3_three_vector::allocate(int min,int max)
00196   {
00197     index_min=min;
00198     index_max=max;
00199     v=new df3_three_variable[max-min+1];
00200     if (v==0)
00201     {
00202       cerr << "error allocating memory in df3_three_vector" << endl;
00203       ad_exit(1);
00204     }
00205     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00206     {
00207       cerr << "Error trying to allocate memory for df3_three_vector"
00208            << endl;;
00209       ad_exit(1);
00210     }
00211     v-=min;
00212   }
00213 
00218   void df3_three_vector::allocate(void)
00219   {
00220     index_min=0;
00221     index_max=-1;
00222     v=0;
00223     shape=0;
00224   }
00225 
00230  dmatrix value(const df3_three_matrix& v)
00231  {
00232    int rmin=v.indexmin();
00233    int rmax=v.indexmax();
00234    dmatrix cm(rmin,rmax);
00235    for (int i=rmin;i<=rmax;i++)
00236    {
00237      int cmin=v(i).indexmin();
00238      int cmax=v(i).indexmax();
00239      cm(i).allocate(cmin,cmax);
00240      for (int j=cmin;j<=cmax;j++)
00241      {
00242        cm(i,j)=value(v(i,j));
00243      }
00244    }
00245    return cm;
00246  }
00247 
00252  df3_three_matrix::df3_three_matrix(const df3_three_matrix& m2)
00253  {
00254    index_min=m2.index_min;
00255    index_max=m2.index_max;
00256    shape=m2.shape;
00257    if (shape)
00258    {
00259      (shape->ncopies)++;
00260    }
00261    v = m2.v;
00262  }
00263 
00268  df3_three_matrix::~df3_three_matrix()
00269  {
00270    deallocate();
00271  }
00272 
00277  void df3_three_matrix::deallocate(void)
00278  {
00279    if (shape)
00280    {
00281      if (shape->ncopies)
00282      {
00283        (shape->ncopies)--;
00284      }
00285      else
00286      {
00287        v = (df3_three_vector*) (shape->get_pointer());
00288        delete [] v;
00289        v=0;
00290        delete shape;
00291        shape=0;
00292      }
00293    }
00294  }
00295 
00300   void df3_three_matrix::initialize(void)
00301   {
00302     int mmin=indexmin();
00303     int mmax=indexmax();
00304     for (int i=mmin;i<=mmax;i++)
00305     {
00306       (*this)(i).initialize();
00307     }
00308   }
00309 
00314   df3_three_matrix::df3_three_matrix(int rmin,int rmax,int cmin,int cmax)
00315   {
00316     index_min=rmin;
00317     index_max=rmax;
00318     v=new df3_three_vector[rmax-rmin+1];
00319     if (v==0)
00320     {
00321       cerr << "error allocating memory in df3_three_matrix" << endl;
00322       ad_exit(1);
00323     }
00324     if ( (shape=new mat_shapex(v)) == NULL)
00325     {
00326       cerr << "Error trying to allocate memory for df3_three_vector"
00327            << endl;;
00328     }
00329     v-=rmin;
00330 
00331     for (int i=rmin;i<=rmax;i++)
00332     {
00333       v[i].allocate(cmin,cmax);
00334     }
00335   }
00336 
00341 df3_three_variable& df3_three_variable::operator-=(const df3_three_variable& v)
00342   {
00343     *get_u() -= *v.get_u();
00344     *get_u_x() -= *v.get_u_x();
00345     *get_u_y() -= *v.get_u_y();
00346     *get_u_z() -= *v.get_u_z();
00347     *get_u_xx() -= *v.get_u_xx();
00348     *get_u_xy() -= *v.get_u_xy();
00349     *get_u_xz() -= *v.get_u_xz();
00350     *get_u_yy() -= *v.get_u_yy();
00351     *get_u_yz() -= *v.get_u_yz();
00352     *get_u_zz() -= *v.get_u_zz();
00353     *get_u_xxx() -= *v.get_u_xxx();
00354     *get_u_xxy() -= *v.get_u_xxy();
00355     *get_u_xxz() -= *v.get_u_xxz();
00356     *get_u_xyy() -= *v.get_u_xyy();
00357     *get_u_xyz() -= *v.get_u_xyz();
00358     *get_u_xzz() -= *v.get_u_xzz();
00359     *get_u_yyy() -= *v.get_u_yyy();
00360     *get_u_yyz() -= *v.get_u_yyz();
00361     *get_u_yzz() -= *v.get_u_yzz();
00362     *get_u_zzz() -= *v.get_u_zzz();
00363     return *this;
00364   }
00365 
00370   df3_three_variable& df3_three_variable::operator -= (double v)
00371   {
00372     *get_u() -= v;
00373     return *this;
00374   }
00375 
00380 df3_three_variable operator-(const df3_three_variable& v)
00381 {
00382   df3_three_variable z;
00383 
00384   *z.get_u() = -(*v.get_u());
00385   *z.get_u_x() = -(*v.get_u_x());
00386   *z.get_u_y() = -(*v.get_u_y());
00387   *z.get_u_z() = -(*v.get_u_z());
00388   *z.get_u_xx() = -(*v.get_u_xx());
00389   *z.get_u_xy() = -(*v.get_u_xy());
00390   *z.get_u_xz() = -(*v.get_u_xz());
00391   *z.get_u_yy() = -(*v.get_u_yy());
00392   *z.get_u_yz() = -(*v.get_u_yz());
00393   *z.get_u_zz() = -(*v.get_u_zz());
00394   *z.get_u_xxx() = -(*v.get_u_xxx());
00395   *z.get_u_xxy() = -(*v.get_u_xxy());
00396   *z.get_u_xxz() = -(*v.get_u_xxz());
00397   *z.get_u_xyy() = -(*v.get_u_xyy());
00398   *z.get_u_xyz() = -(*v.get_u_xyz());
00399   *z.get_u_xzz() = -(*v.get_u_xzz());
00400   *z.get_u_yyy() = -(*v.get_u_yyy());
00401   *z.get_u_yyz() = -(*v.get_u_yyz());
00402   *z.get_u_yzz() = -(*v.get_u_yzz());
00403   *z.get_u_zzz() = -(*v.get_u_zzz());
00404 
00405   return z;
00406 }
00407 
00412 df3_three_variable& df3_three_variable::operator+=(const df3_three_variable& v)
00413   {
00414     *get_u() += *v.get_u();
00415     *get_u_x() += *v.get_u_x();
00416     *get_u_y() += *v.get_u_y();
00417     *get_u_z() += *v.get_u_z();
00418     *get_u_xx() += *v.get_u_xx();
00419     *get_u_xy() += *v.get_u_xy();
00420     *get_u_xz() += *v.get_u_xz();
00421     *get_u_yy() += *v.get_u_yy();
00422     *get_u_yz() += *v.get_u_yz();
00423     *get_u_zz() += *v.get_u_zz();
00424     *get_u_xxx() += *v.get_u_xxx();
00425     *get_u_xxy() += *v.get_u_xxy();
00426     *get_u_xxz() += *v.get_u_xxz();
00427     *get_u_xyy() += *v.get_u_xyy();
00428     *get_u_xyz() += *v.get_u_xyz();
00429     *get_u_xzz() += *v.get_u_xzz();
00430     *get_u_yyy() += *v.get_u_yyy();
00431     *get_u_yyz() += *v.get_u_yyz();
00432     *get_u_yzz() += *v.get_u_yzz();
00433     *get_u_zzz() += *v.get_u_zzz();
00434     return *this;
00435   }
00436 
00441   df3_three_variable& df3_three_variable::operator *= (double v)
00442   {
00443     *get_u() *= v;
00444     *get_u_x() *= v;
00445     *get_u_y() *= v;
00446     *get_u_z() *= v;
00447     *get_u_xx() *= v;
00448     *get_u_xy() *= v;
00449     *get_u_xz() *= v;
00450     *get_u_yy() *= v;
00451     *get_u_yz() *= v;
00452     *get_u_zz() *= v;
00453     *get_u_xxx() *= v;
00454     *get_u_xxy() *= v;
00455     *get_u_xxz() *= v;
00456     *get_u_xyy() *= v;
00457     *get_u_xyz() *= v;
00458     *get_u_xzz() *= v;
00459     *get_u_yyy() *= v;
00460     *get_u_yyz() *= v;
00461     *get_u_yzz() *= v;
00462     *get_u_zzz() *= v;
00463     return *this;
00464   }
00465 
00470 df3_three_variable& df3_three_variable::operator*=(const df3_three_variable& v)
00471   {
00472     df3_three_variable x=*this * v;
00473     *this=x;
00474     return *this;
00475   }
00476 
00481 df3_three_variable& df3_three_variable::operator/=(const df3_three_variable& v)
00482   {
00483     df3_three_variable x=*this / v;
00484     *this=x;
00485     return *this;
00486   }
00487 
00492   df3_three_variable& df3_three_variable::operator += (double v)
00493   {
00494     *get_u() += v;
00495 
00496     return *this;
00497   }
00498 
00503 void set_derivatives(df3_three_variable& z,const df3_three_variable& x,
00504   double u, double zp,double zp2,double zp3)
00505 {
00506     //*z.get_u() = u;
00507 
00508     *z.get_u_x() = zp* *x.get_u_x();
00509 
00510     *z.get_u_y() = zp* *x.get_u_y();
00511 
00512     *z.get_u_z() = zp* *x.get_u_z();
00513 
00514     *z.get_u_xx() = zp2 * square(*x.get_u_x())
00515                    + zp * *x.get_u_xx();
00516 
00517     *z.get_u_xy() = zp2 * *x.get_u_x() * *x.get_u_y()
00518                    + zp * *x.get_u_xy();
00519 
00520     *z.get_u_xz() = zp2 * *x.get_u_x() * *x.get_u_z()
00521                    + zp * *x.get_u_xz();
00522 
00523     *z.get_u_yy() = zp2 * square(*x.get_u_y())
00524                    + zp * *x.get_u_yy();
00525 
00526     *z.get_u_yz() = zp2 * *x.get_u_y()* *x.get_u_z()
00527                    + zp * *x.get_u_yz();
00528 
00529     *z.get_u_zz() = zp2 * square(*x.get_u_z())
00530                    + zp * *x.get_u_zz();
00531 
00532     *z.get_u_xxx() = zp3 * cube(*x.get_u_x())
00533                    + 3.0 * zp2 * *x.get_u_x() * *x.get_u_xx()
00534                    + zp * *x.get_u_xxx();
00535 
00536     *z.get_u_xxy() = zp3 * square(*x.get_u_x())* *x.get_u_y()
00537                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xy()
00538                    + zp2 * *x.get_u_y() * *x.get_u_xx()
00539                    + zp * *x.get_u_xxy();
00540 
00541     *z.get_u_xxz() = zp3 * square(*x.get_u_x())* *x.get_u_z()
00542                    + 2.0* zp2 * *x.get_u_x() * *x.get_u_xz()
00543                    + zp2 * *x.get_u_z() * *x.get_u_xx()
00544                    + zp * *x.get_u_xxz();
00545 
00546     *z.get_u_xyy() = zp3 * square(*x.get_u_y())* *x.get_u_x()
00547                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_xy()
00548                    + zp2 * *x.get_u_x() * *x.get_u_yy()
00549                    + zp * *x.get_u_xyy();
00550 
00551     *z.get_u_xyz() = zp3 * *x.get_u_y()* *x.get_u_x()* *x.get_u_z()
00552                    + zp2 * *x.get_u_z() * *x.get_u_xy()
00553                    + zp2 * *x.get_u_y() * *x.get_u_xz()
00554                    + zp2 * *x.get_u_x() * *x.get_u_yz()
00555                    + zp * *x.get_u_xyz();
00556 
00557     *z.get_u_xzz() = zp3 * square(*x.get_u_z())* *x.get_u_x()
00558                    + 2.0* zp2 * *x.get_u_z() * *x.get_u_xz()
00559                    + zp2 * *x.get_u_x() * *x.get_u_zz()
00560                    + zp * *x.get_u_xzz();
00561 
00562 
00563     *z.get_u_yyy() = zp3 * cube(*x.get_u_y())
00564                    + 3.0 * zp2 * *x.get_u_y() * *x.get_u_yy()
00565                    + zp * *x.get_u_yyy();
00566     *z.get_u_yyz() = zp3 * square(*x.get_u_y())* *x.get_u_z()
00567                    + 2.0* zp2 * *x.get_u_y() * *x.get_u_yz()
00568                    + zp2 * *x.get_u_z() * *x.get_u_yy()
00569                    + zp * *x.get_u_yyz();
00570 
00571     *z.get_u_yzz() = zp3 * square(*x.get_u_z())* *x.get_u_y()
00572                    + 2.0* zp2 * *x.get_u_z() * *x.get_u_yz()
00573                    + zp2 * *x.get_u_y() * *x.get_u_zz()
00574                    + zp * *x.get_u_yzz();
00575 
00576     *z.get_u_zzz() = zp3 * cube(*x.get_u_z())
00577                    + 3.0 * zp2 * *x.get_u_z() * *x.get_u_zz()
00578                    + zp * *x.get_u_zzz();
00579 }
00580 
00585 void set_derivatives( df3_three_variable& z, const df3_three_variable& x,
00586   const df3_three_variable& y, double u,
00587   double f_u,double f_v,double f_uu,double f_uv,double f_vv,double f_uuu,
00588   double f_uuv,double f_uvv,double f_vvv)
00589 {
00590     *z.get_u() = u;
00591 
00592     *z.get_u_x() = f_u* *x.get_u_x()
00593                  + f_v* *y.get_u_x();
00594 
00595     *z.get_u_y() = f_u* *x.get_u_y()
00596                  + f_v* *y.get_u_y();
00597 
00598     *z.get_u_z() = f_u* *x.get_u_z()
00599                  + f_v* *y.get_u_z();
00600 
00601     *z.get_u_xx() = f_uu * square(*x.get_u_x())
00602                   + f_u  * *x.get_u_xx()
00603                   + f_vv * square(*y.get_u_x())
00604                   + f_v  * *y.get_u_xx()
00605             + 2.0 * f_uv * *x.get_u_x() * *y.get_u_x();
00606 
00607     *z.get_u_xy() = f_uu * *x.get_u_x() * *x.get_u_y()
00608                   + f_u  * *x.get_u_xy()
00609                   + f_vv *  *y.get_u_x() * *y.get_u_y()
00610                   + f_v  * *y.get_u_xy()
00611                   + f_uv * (*x.get_u_x() * *y.get_u_y()
00612                          +  *x.get_u_y() * *y.get_u_x());
00613 
00614     *z.get_u_xz() = f_uu * *x.get_u_x() * *x.get_u_z()
00615                   + f_u  * *x.get_u_xz()
00616                   + f_vv *  *y.get_u_x() * *y.get_u_z()
00617                   + f_v  * *y.get_u_xz()
00618                   + f_uv * (*x.get_u_x() * *y.get_u_z()
00619                          +  *x.get_u_z() * *y.get_u_x());
00620 
00621     *z.get_u_yy() = f_uu * square(*x.get_u_y())
00622                   + f_u  * *x.get_u_yy()
00623                   + f_vv * square(*y.get_u_y())
00624                   + f_v  * *y.get_u_yy()
00625             + 2.0 * f_uv * *x.get_u_y() * *y.get_u_y();
00626 
00627     *z.get_u_yz() = f_uu * *x.get_u_y() * *x.get_u_z()
00628                   + f_u  * *x.get_u_yz()
00629                   + f_vv *  *y.get_u_y() * *y.get_u_z()
00630                   + f_v  * *y.get_u_yz()
00631                   + f_uv * (*x.get_u_y() * *y.get_u_z()
00632                          +  *x.get_u_z() * *y.get_u_y());
00633 
00634     *z.get_u_zz() = f_uu * *x.get_u_z() * *x.get_u_z()
00635                   + f_u  * *x.get_u_zz()
00636                   + f_vv *  *y.get_u_z() * *y.get_u_z()
00637                   + f_v  * *y.get_u_zz()
00638                   + f_uv * (*x.get_u_z() * *y.get_u_z()
00639                          +  *x.get_u_z() * *y.get_u_z());
00640 
00641 
00642     *z.get_u_xxx() = f_uuu * cube(*x.get_u_x())
00643                    + 3.0 * f_uu * *x.get_u_x() * *x.get_u_xx()
00644                    + f_u * *x.get_u_xxx()
00645                    + f_vvv * cube(*y.get_u_x())
00646                    + 3.0 * f_vv * *y.get_u_x() * *y.get_u_xx()
00647                    + f_v * *y.get_u_xxx()
00648                    + 3.0 * f_uuv * square(*x.get_u_x()) * *y.get_u_x()
00649                    + 3.0 * f_uvv * *x.get_u_x()* square(*y.get_u_x())
00650                    + 3.0* f_uv * *x.get_u_xx() * *y.get_u_x()
00651                    + 3.0* f_uv * *x.get_u_x() * *y.get_u_xx();
00652 
00653     *z.get_u_xxy() = f_uuu * square(*x.get_u_x())* *x.get_u_y()
00654                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xy()
00655                    + f_uu * *x.get_u_y() * *x.get_u_xx()
00656                    + f_u * *x.get_u_xxy()
00657                    + f_vvv * square(*y.get_u_x())* *y.get_u_y()
00658                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xy()
00659                    + f_vv * *y.get_u_xx() * *y.get_u_y()
00660                    + f_v * *y.get_u_xxy()
00661                    + f_uuv * square(*x.get_u_x()) * *y.get_u_y()
00662                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_y() * *y.get_u_x()
00663                    + f_uvv * *x.get_u_y()* square(*y.get_u_x())
00664                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_y()
00665                    + f_uv * *x.get_u_xx() * *y.get_u_y()
00666                    + f_uv * *x.get_u_y() * *y.get_u_xx()
00667                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_x()
00668                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xy();
00669 
00670     *z.get_u_xxz() = f_uuu * square(*x.get_u_x())* *x.get_u_z()
00671                    + 2.0 * f_uu * *x.get_u_x() * *x.get_u_xz()
00672                    + f_uu * *x.get_u_z() * *x.get_u_xx()
00673                    + f_u * *x.get_u_xxz()
00674                    + f_vvv * square(*y.get_u_x())* *y.get_u_z()
00675                    + 2.0 * f_vv * *y.get_u_x() * *y.get_u_xz()
00676                    + f_vv * *y.get_u_xx() * *y.get_u_z()
00677                    + f_v * *y.get_u_xxz()
00678                    + f_uuv * square(*x.get_u_x()) * *y.get_u_z()
00679                    + 2.0* f_uuv * *x.get_u_x() * *x.get_u_z() * *y.get_u_x()
00680                    + f_uvv * *x.get_u_z()* square(*y.get_u_x())
00681                    + 2.0 * f_uvv * *x.get_u_x()* *y.get_u_x() * *y.get_u_z()
00682                    + f_uv * *x.get_u_xx() * *y.get_u_z()
00683                    + f_uv * *x.get_u_z() * *y.get_u_xx()
00684                    + 2.0* f_uv * *x.get_u_xz() * *y.get_u_x()
00685                    + 2.0* f_uv * *x.get_u_x() * *y.get_u_xz();
00686 
00687     *z.get_u_xyy() = f_uuu * square(*x.get_u_y())* *x.get_u_x()
00688                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_xy()
00689                    + f_uu * *x.get_u_x() * *x.get_u_yy()
00690                    + f_u * *x.get_u_xyy()
00691                    + f_vvv * square(*y.get_u_y())* *y.get_u_x()
00692                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_xy()
00693                    + f_vv * *y.get_u_yy() * *y.get_u_x()
00694                    + f_v * *y.get_u_xyy()
00695                    + f_uuv * square(*x.get_u_y()) * *y.get_u_x()
00696                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_x() * *y.get_u_y()
00697                    + f_uvv * *x.get_u_x()* square(*y.get_u_y())
00698                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_x()
00699                    + f_uv * *x.get_u_yy() * *y.get_u_x()
00700                    + f_uv * *x.get_u_x() * *y.get_u_yy()
00701                    + 2.0* f_uv * *x.get_u_xy() * *y.get_u_y()
00702                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_xy();
00703 
00704     *z.get_u_xyz() = f_uuu * *x.get_u_y()* *x.get_u_x()* *x.get_u_z()
00705                    + f_uuv * *x.get_u_y() * *x.get_u_x()* *y.get_u_z()
00706                    + f_uu * *x.get_u_y() * *x.get_u_xz()
00707                    + f_uu * *x.get_u_x() * *x.get_u_yz()
00708 
00709                    + f_uuv * *y.get_u_y() * *x.get_u_x()* *x.get_u_z()
00710                    + f_uvv * *y.get_u_y() * *x.get_u_x()* *y.get_u_z()
00711                    + f_uv * *x.get_u_xz() * *y.get_u_y()
00712                    + f_uv * *y.get_u_yz() * *x.get_u_x()
00713 
00714                    + f_uu * *x.get_u_xy() * *x.get_u_z()
00715                    + f_uv * *x.get_u_xy() * *y.get_u_z()
00716                    + f_u * *x.get_u_xyz()
00717 
00718                    + f_uuv * *x.get_u_y() * *y.get_u_x()* *x.get_u_z()
00719                    + f_uvv * *x.get_u_y() * *y.get_u_x()* *y.get_u_z()
00720                    + f_uv * *y.get_u_x() * *x.get_u_yz()
00721                    + f_uv * *x.get_u_y() * *y.get_u_xz()
00722 
00723                    + f_uvv * *y.get_u_y() * *y.get_u_x()* *x.get_u_z()
00724                    + f_vvv * *y.get_u_y() * *y.get_u_x()* *y.get_u_z()
00725 
00726                    + f_vv * *y.get_u_xz() * *y.get_u_y()
00727                    + f_vv * *y.get_u_yz() * *y.get_u_x()
00728                    + f_vv * *y.get_u_xy() * *y.get_u_z()
00729                    + f_uv * *y.get_u_xy() * *x.get_u_z()
00730                    + f_v * *y.get_u_xyz();
00731 
00732     *z.get_u_xzz() = f_uuu * square(*x.get_u_z())* *x.get_u_x()
00733                    + 2.0 * f_uu * *x.get_u_z() * *x.get_u_xz()
00734                    + f_uu * *x.get_u_x() * *x.get_u_zz()
00735                    + f_u * *x.get_u_xzz()
00736                    + f_vvv * square(*y.get_u_z())* *y.get_u_x()
00737                    + 2.0 * f_vv * *y.get_u_z() * *y.get_u_xz()
00738                    + f_vv * *y.get_u_zz() * *y.get_u_x()
00739                    + f_v * *y.get_u_xzz()
00740                    + f_uuv * square(*x.get_u_z()) * *y.get_u_x()
00741                    + 2.0* f_uuv * *x.get_u_z() * *x.get_u_x() * *y.get_u_z()
00742                    + f_uvv * *x.get_u_x()* square(*y.get_u_z())
00743                    + 2.0 * f_uvv * *x.get_u_z()* *y.get_u_z() * *y.get_u_x()
00744                    + f_uv * *x.get_u_zz() * *y.get_u_x()
00745                    + f_uv * *x.get_u_x() * *y.get_u_zz()
00746                    + 2.0* f_uv * *x.get_u_xz() * *y.get_u_z()
00747                    + 2.0* f_uv * *x.get_u_z() * *y.get_u_xz();
00748 
00749     *z.get_u_yyy() = f_uuu * cube(*x.get_u_y())
00750                    + 3.0 * f_uu * *x.get_u_y() * *x.get_u_yy()
00751                    + f_u * *x.get_u_yyy()
00752                    + f_vvv * cube(*y.get_u_y())
00753                    + 3.0 * f_vv * *y.get_u_y() * *y.get_u_yy()
00754                    + f_v * *y.get_u_yyy()
00755                    + 3.0 * f_uuv * square(*x.get_u_y()) * *y.get_u_y()
00756                    + 3.0 * f_uvv * *x.get_u_y()* square(*y.get_u_y())
00757                    + 3.0 * f_uv * *x.get_u_yy() * *y.get_u_y()
00758                    + 3.0 * f_uv * *x.get_u_y() * *y.get_u_yy();
00759 
00760     *z.get_u_yyz() = f_uuu * square(*x.get_u_y())* *x.get_u_z()
00761                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_yz()
00762                    + f_uu * *x.get_u_z() * *x.get_u_yy()
00763                    + f_u * *x.get_u_yyz()
00764                    + f_vvv * square(*y.get_u_y())* *y.get_u_z()
00765                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_yz()
00766                    + f_vv * *y.get_u_yy() * *y.get_u_z()
00767                    + f_v * *y.get_u_yyz()
00768                    + f_uuv * square(*x.get_u_y()) * *y.get_u_z()
00769                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_z() * *y.get_u_y()
00770                    + f_uvv * *x.get_u_z()* square(*y.get_u_y())
00771                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_z()
00772                    + f_uv * *x.get_u_yy() * *y.get_u_z()
00773                    + f_uv * *x.get_u_z() * *y.get_u_yy()
00774                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_y()
00775                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_yz();
00776 
00777     *z.get_u_yyz() = f_uuu * square(*x.get_u_y())* *x.get_u_z()
00778                    + 2.0 * f_uu * *x.get_u_y() * *x.get_u_yz()
00779                    + f_uu * *x.get_u_z() * *x.get_u_yy()
00780                    + f_u * *x.get_u_yyz()
00781                    + f_vvv * square(*y.get_u_y())* *y.get_u_z()
00782                    + 2.0 * f_vv * *y.get_u_y() * *y.get_u_yz()
00783                    + f_vv * *y.get_u_yy() * *y.get_u_z()
00784                    + f_v * *y.get_u_yyz()
00785                    + f_uuv * square(*x.get_u_y()) * *y.get_u_z()
00786                    + 2.0* f_uuv * *x.get_u_y() * *x.get_u_z() * *y.get_u_y()
00787                    + f_uvv * *x.get_u_z()* square(*y.get_u_y())
00788                    + 2.0 * f_uvv * *x.get_u_y()* *y.get_u_y() * *y.get_u_z()
00789                    + f_uv * *x.get_u_yy() * *y.get_u_z()
00790                    + f_uv * *x.get_u_z() * *y.get_u_yy()
00791                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_y()
00792                    + 2.0* f_uv * *x.get_u_y() * *y.get_u_yz();
00793 
00794     *z.get_u_yzz() = f_uuu * square(*x.get_u_z())* *x.get_u_y()
00795                    + 2.0 * f_uu * *x.get_u_z() * *x.get_u_yz()
00796                    + f_uu * *x.get_u_y() * *x.get_u_zz()
00797                    + f_u * *x.get_u_yzz()
00798                    + f_vvv * square(*y.get_u_z())* *y.get_u_y()
00799                    + 2.0 * f_vv * *y.get_u_z() * *y.get_u_yz()
00800                    + f_vv * *y.get_u_zz() * *y.get_u_y()
00801                    + f_v * *y.get_u_yzz()
00802                    + f_uuv * square(*x.get_u_z()) * *y.get_u_y()
00803                    + 2.0* f_uuv * *x.get_u_z() * *x.get_u_y() * *y.get_u_z()
00804                    + f_uvv * *x.get_u_y()* square(*y.get_u_z())
00805                    + 2.0 * f_uvv * *x.get_u_z()* *y.get_u_z() * *y.get_u_y()
00806                    + f_uv * *x.get_u_zz() * *y.get_u_y()
00807                    + f_uv * *x.get_u_y() * *y.get_u_zz()
00808                    + 2.0* f_uv * *x.get_u_yz() * *y.get_u_z()
00809                    + 2.0* f_uv * *x.get_u_z() * *y.get_u_yz();
00810 
00811     *z.get_u_zzz() = f_uuu * cube(*x.get_u_z())
00812                    + 3.0 * f_uu * *x.get_u_z() * *x.get_u_zz()
00813                    + f_u * *x.get_u_zzz()
00814                    + f_vvv * cube(*y.get_u_z())
00815                    + 3.0 * f_vv * *y.get_u_z() * *y.get_u_zz()
00816                    + f_v * *y.get_u_zzz()
00817                    + 3.0 * f_uuv * square(*x.get_u_z()) * *y.get_u_z()
00818                    + 3.0 * f_uvv * *x.get_u_z()* square(*y.get_u_z())
00819                    + 3.0 * f_uv * *x.get_u_zz() * *y.get_u_z()
00820                    + 3.0 * f_uv * *x.get_u_z() * *y.get_u_zz();
00821 }
00822 
00827   df3_three_variable sqrt(const df3_three_variable& x)
00828   {
00829     df3_three_variable z;
00830     double u=::sqrt(*x.get_u());
00831     *z.get_u()=u;
00832     double xinv=1.0/(*x.get_u());
00833     double zp=0.5/u;
00834     double zp2=-0.5*zp*xinv;
00835     double zp3=-1.5*zp2*xinv;
00836 
00837 
00838     set_derivatives(z,x,u,zp,zp2,zp3);
00839 
00840     return z;
00841   }
00842 
00847   df3_three_variable atan(const df3_three_variable& x)
00848   {
00849     df3_three_variable z;
00850     double cx=value(x);
00851     double d=1.0/(1+square(cx));
00852     double d2=square(d);
00853     double u=::atan(cx);
00854     *z.get_u()=u;
00855     double zp=d;
00856     double zp2=-2.0*cx*d2;
00857     double zp3=-2.0*d2+8*cx*cx*d*d2;
00858 
00859     set_derivatives(z,x,u,zp,zp2,zp3);
00860     return z;
00861   }
00862 
00867   df3_three_variable pow(const df3_three_variable& x,
00868     const df3_three_variable& y)
00869   {
00870     df3_three_variable z;
00871     z=exp(y*log(x));
00872     return z;
00873   }
00874 
00879   df3_three_variable square(const df3_three_variable& x)
00880   {
00881     df3_three_variable z;
00882     double u=value(x);
00883     *z.get_u()=u*u;
00884     double zp=2.0*u;
00885     double zp2=2.0;
00886     double zp3=0.0;
00887 
00888     set_derivatives(z,x,u,zp,zp2,zp3);
00889     return z;
00890   }
00891 
00896   df3_three_variable tan(const df3_three_variable& x)
00897   {
00898     df3_three_variable z;
00899     double u=::tan(*x.get_u());
00900     *z.get_u()=u;
00901     double v=1.0/::cos(*x.get_u());
00902     double w=::sin(*x.get_u());
00903     double v2=v*v;
00904     double zp=v2;
00905     double zp2=2.0*w*v2*v;
00906     double zp3=(4.0*w*w+2.0)*v2*v2;
00907 
00908     set_derivatives(z,x,u,zp,zp2,zp3);
00909     return z;
00910   }
00911 
00916   df3_three_variable sin(const df3_three_variable& x)
00917   {
00918     df3_three_variable z;
00919     double u=::sin(*x.get_u());
00920     *z.get_u()=u;
00921     double zp=::cos(*x.get_u());
00922     double zp2=-u;
00923     double zp3=-zp;
00924 
00925     set_derivatives(z,x,u,zp,zp2,zp3);
00926     return z;
00927   }
00928 
00933   df3_three_variable log(const df3_three_variable& x)
00934   {
00935     df3_three_variable z;
00936     double u=::log(*x.get_u());
00937     *z.get_u()=u;
00938     double zp=1/(*x.get_u());
00939     double zp2=-zp*zp;
00940     double zp3=-2.0*zp*zp2;
00941 
00942     set_derivatives(z,x,u,zp,zp2,zp3);
00943     return z;
00944   }
00945 
00950   df3_three_variable exp(const df3_three_variable& x)
00951   {
00952     df3_three_variable z;
00953     double u=::exp(*x.get_u());
00954     *z.get_u()=u;
00955     double zp=u;
00956     double zp2=u;
00957     double zp3=u;
00958 
00959     set_derivatives(z,x,u,zp,zp2,zp3);
00960     return z;
00961   }
00962 
00967   df3_three_variable inv(const df3_three_variable& x)
00968   {
00969     df3_three_variable z;
00970     double xinv=1.0/(*x.get_u());
00971     *z.get_u()=xinv;
00972     double zp=-xinv*xinv;
00973     double zp2=-2.0*zp*xinv;
00974     double zp3=-3.0*zp2*xinv;
00975 
00976     set_derivatives(z,x,xinv,zp,zp2,zp3);
00977 
00978     return z;
00979   }
00980 
00985 df3_three_variable& df3_three_variable::operator=(const df3_three_variable& v)
00986   {
00987     *get_u() =  *v.get_u();
00988     *get_u_x() = *v.get_u_x();
00989     *get_u_y() = *v.get_u_y();
00990     *get_u_z() = *v.get_u_z();
00991     *get_u_xx() = *v.get_u_xx();
00992     *get_u_xy() = *v.get_u_xy();
00993     *get_u_xz() = *v.get_u_xz();
00994     *get_u_yy() = *v.get_u_yy();
00995     *get_u_yz() = *v.get_u_yz();
00996     *get_u_zz() = *v.get_u_zz();
00997     *get_u_xxx() = *v.get_u_xxx();
00998     *get_u_xxy() = *v.get_u_xxy();
00999     *get_u_xxz() = *v.get_u_xxz();
01000     *get_u_xyy() = *v.get_u_xyy();
01001     *get_u_xyz() = *v.get_u_xyz();
01002     *get_u_xzz() = *v.get_u_xzz();
01003     *get_u_yyy() = *v.get_u_yyy();
01004     *get_u_yyz() = *v.get_u_yyz();
01005     *get_u_yzz() = *v.get_u_yzz();
01006     *get_u_zzz() = *v.get_u_zzz();
01007 
01008 
01009 
01010 
01011     return *this;
01012   }
01013 
01018   df3_three_variable& df3_three_variable::operator = (double x)
01019   {
01020     *get_u() =x;
01021     *get_u_x() =0.0;
01022     *get_u_y() =0.0;
01023     *get_u_z() =0.0;
01024     *get_u_xx() =0.0;
01025     *get_u_xy() =0.0;
01026     *get_u_xz() =0.0;
01027     *get_u_yy() =0.0;
01028     *get_u_yz() =0.0;
01029     *get_u_zz() =0.0;
01030     *get_u_xxx() =0.0;
01031     *get_u_xxy() =0.0;
01032     *get_u_xxz() =0.0;
01033     *get_u_xyy() =0.0;
01034     *get_u_xyz() =0.0;
01035     *get_u_xzz() =0.0;
01036     *get_u_yyy() =0.0;
01037     *get_u_yyz() =0.0;
01038     *get_u_yzz() =0.0;
01039     *get_u_zzz() =0.0;
01040 
01041     return *this;
01042   }
01043 
01048   df3_three_variable operator * (const df3_three_variable& x,
01049     const df3_three_variable& y)
01050   {
01051     df3_three_variable z;
01052     double u= *x.get_u() * *y.get_u();
01053     *z.get_u() = u;
01054     double f_u=*y.get_u();
01055     double f_v=*x.get_u();
01056     double f_uu=0.0;
01057     double f_uv=1.0;
01058     double f_vv=0.0;
01059     double f_uuu=0.0;
01060     double f_uuv=0.0;
01061     double f_uvv=0.0;
01062     double f_vvv=0.0;
01063     set_derivatives(z,x,y,u,
01064       f_u, f_v,
01065       f_uu, f_uv, f_vv,
01066       f_uuu, f_uuv, f_uvv, f_vvv);
01067     return z;
01068   }
01069 
01074   df3_three_variable operator * (double x,
01075     const df3_three_variable& v)
01076   {
01077     df3_three_variable z;
01078 
01079     *z.get_u() = x * *v.get_u();
01080     *z.get_u_x() = x * *v.get_u_x();
01081     *z.get_u_y() = x * *v.get_u_y();
01082     *z.get_u_z() = x * *v.get_u_z();
01083     *z.get_u_xx() = x * *v.get_u_xx();
01084     *z.get_u_xy() = x * *v.get_u_xy();
01085     *z.get_u_xz() = x * *v.get_u_xz();
01086     *z.get_u_yy() = x * *v.get_u_yy();
01087     *z.get_u_yz() = x * *v.get_u_yz();
01088     *z.get_u_zz() = x * *v.get_u_zz();
01089     *z.get_u_xxx() = x * *v.get_u_xxx();
01090     *z.get_u_xxy() = x * *v.get_u_xxy();
01091     *z.get_u_xxz() = x * *v.get_u_xxz();
01092     *z.get_u_xyy() = x * *v.get_u_xyy();
01093     *z.get_u_xyz() = x * *v.get_u_xyz();
01094     *z.get_u_xzz() = x * *v.get_u_xzz();
01095     *z.get_u_yyy() = x * *v.get_u_yyy();
01096     *z.get_u_yyz() = x * *v.get_u_yyz();
01097     *z.get_u_yzz() = x * *v.get_u_yzz();
01098     *z.get_u_zzz() = x * *v.get_u_zzz();
01099 
01100     return z;
01101   }
01102 
01107   df3_three_variable fabs(const df3_three_variable& v)
01108   {
01109     df3_three_variable z;
01110     if (value(v)>=0)
01111     {
01112       *z.get_u() = *v.get_u();
01113       *z.get_u_x() = *v.get_u_x();
01114       *z.get_u_y() = *v.get_u_y();
01115       *z.get_u_z() = *v.get_u_z();
01116       *z.get_u_xx() = *v.get_u_xx();
01117       *z.get_u_xy() = *v.get_u_xy();
01118       *z.get_u_xz() = *v.get_u_xz();
01119       *z.get_u_yy() = *v.get_u_yy();
01120       *z.get_u_yz() = *v.get_u_yz();
01121       *z.get_u_zz() = *v.get_u_zz();
01122       *z.get_u_xxx() = *v.get_u_xxx();
01123       *z.get_u_xxy() = *v.get_u_xxy();
01124       *z.get_u_xxz() = *v.get_u_xxz();
01125       *z.get_u_xyy() = *v.get_u_xyy();
01126       *z.get_u_xyz() = *v.get_u_xyz();
01127       *z.get_u_xzz() = *v.get_u_xzz();
01128       *z.get_u_yyy() = *v.get_u_yyy();
01129       *z.get_u_yyz() = *v.get_u_yyz();
01130       *z.get_u_yzz() = *v.get_u_yzz();
01131       *z.get_u_zzz() = *v.get_u_zzz();
01132     }
01133     else
01134     {
01135       *z.get_u() = -*v.get_u();
01136       *z.get_u_x() = -*v.get_u_x();
01137       *z.get_u_y() = -*v.get_u_y();
01138       *z.get_u_z() = -*v.get_u_z();
01139       *z.get_u_xx() = -*v.get_u_xx();
01140       *z.get_u_xy() = -*v.get_u_xy();
01141       *z.get_u_xz() = -*v.get_u_xz();
01142       *z.get_u_yy() = -*v.get_u_yy();
01143       *z.get_u_yz() = -*v.get_u_yz();
01144       *z.get_u_zz() = -*v.get_u_zz();
01145       *z.get_u_xxx() = -*v.get_u_xxx();
01146       *z.get_u_xxy() = -*v.get_u_xxy();
01147       *z.get_u_xxz() = -*v.get_u_xxz();
01148       *z.get_u_xyy() = -*v.get_u_xyy();
01149       *z.get_u_xyz() = -*v.get_u_xyz();
01150       *z.get_u_xzz() = -*v.get_u_xzz();
01151       *z.get_u_yyy() = -*v.get_u_yyy();
01152       *z.get_u_yyz() = -*v.get_u_yyz();
01153       *z.get_u_yzz() = -*v.get_u_yzz();
01154       *z.get_u_zzz() = -*v.get_u_zzz();
01155     }
01156 
01157     return z;
01158   }
01159 
01164   df3_three_variable operator * (const df3_three_variable& v,
01165     double x)
01166   {
01167     df3_three_variable z;
01168 
01169     *z.get_u() = x * *v.get_u();
01170     *z.get_u_x() = x * *v.get_u_x();
01171     *z.get_u_y() = x * *v.get_u_y();
01172     *z.get_u_z() = x * *v.get_u_z();
01173     *z.get_u_xx() = x * *v.get_u_xx();
01174     *z.get_u_xy() = x * *v.get_u_xy();
01175     *z.get_u_xz() = x * *v.get_u_xz();
01176     *z.get_u_yy() = x * *v.get_u_yy();
01177     *z.get_u_yz() = x * *v.get_u_yz();
01178     *z.get_u_zz() = x * *v.get_u_zz();
01179     *z.get_u_xxx() = x * *v.get_u_xxx();
01180     *z.get_u_xxy() = x * *v.get_u_xxy();
01181     *z.get_u_xxz() = x * *v.get_u_xxz();
01182     *z.get_u_xyy() = x * *v.get_u_xyy();
01183     *z.get_u_xyz() = x * *v.get_u_xyz();
01184     *z.get_u_xzz() = x * *v.get_u_xzz();
01185     *z.get_u_yyy() = x * *v.get_u_yyy();
01186     *z.get_u_yyz() = x * *v.get_u_yyz();
01187     *z.get_u_yzz() = x * *v.get_u_yzz();
01188     *z.get_u_zzz() = x * *v.get_u_zzz();
01189 
01190     return z;
01191   }
01192 
01197   df3_three_variable operator / (const df3_three_variable& x,
01198     double y)
01199   {
01200     double u=1/y;
01201     return x*u;
01202   }
01203 
01208   df3_three_variable operator / (const df3_three_variable& x,
01209     const df3_three_variable& y)
01210   {
01211     df3_three_variable u=inv(y);
01212     return x*u;
01213   }
01214 
01219   df3_three_variable operator / (const double x,
01220     const df3_three_variable& y)
01221   {
01222     df3_three_variable u=inv(y);
01223     return x*u;
01224   }
01225 
01226  /*
01227   df3_three_variable operator / (const df3_three_variable& x,
01228     const df3_three_variable& y)
01229   {
01230     df3_three_variable z;
01231     double yinv =  1.0 / (*y.get_u());
01232     double yinv2 = yinv * yinv;
01233     double yinv3 = yinv * yinv2;
01234     doubl yd = *y.get_udot();
01235 
01236     // *z.get_u() = *x.get_u() /  *y.get_u();
01237     *z.get_u() = *x.get_u() * yinv;
01238 
01239     *z.get_udot() =  - (*x.get_u()) * yinv2 * yd
01240                   + *x.get_udot() * yinv;
01241 
01242     *z.get_udot2() = *x.get_udot2() * yinv
01243                    - 2.0 * *x.get_udot() * yd * yinv2
01244                    + 2.0 * *x.get_u() * yinv3  * yd *yd
01245                    -  *x.get_u() * yinv2 * y.get_udot2();
01246 
01247     *z.get_udot3() = *x.get_udot3() * yinv
01248                    + 3.0 * *x.get_udot2() * *y.get_udot()
01249                    + 3.0 * *x.get_udot() * *y.get_udot2()
01250                    +  *x.get_u() * *y.get_udot3();
01251   }
01252  */
01253 
01258   df3_three_variable operator + (const double x,const df3_three_variable& v)
01259   {
01260     df3_three_variable z;
01261     *z.get_u() = x + *v.get_u();
01262     *z.get_u_x() = *v.get_u_x();
01263     *z.get_u_y() = *v.get_u_y();
01264     *z.get_u_z() = *v.get_u_z();
01265     *z.get_u_xx() = *v.get_u_xx();
01266     *z.get_u_xy() = *v.get_u_xy();
01267     *z.get_u_xz() = *v.get_u_xz();
01268     *z.get_u_yy() = *v.get_u_yy();
01269     *z.get_u_yz() = *v.get_u_yz();
01270     *z.get_u_zz() = *v.get_u_zz();
01271     *z.get_u_xxx() = *v.get_u_xxx();
01272     *z.get_u_xxy() = *v.get_u_xxy();
01273     *z.get_u_xxz() = *v.get_u_xxz();
01274     *z.get_u_xyy() = *v.get_u_xyy();
01275     *z.get_u_xyz() = *v.get_u_xyz();
01276     *z.get_u_xzz() = *v.get_u_xzz();
01277     *z.get_u_yyy() = *v.get_u_yyy();
01278     *z.get_u_yyz() = *v.get_u_yyz();
01279     *z.get_u_yzz() = *v.get_u_yzz();
01280     *z.get_u_zzz() = *v.get_u_zzz();
01281 
01282     return z;
01283   }
01284 
01289   df3_three_variable operator + (const df3_three_variable& v,const double x)
01290   {
01291     df3_three_variable z;
01292 
01293     *z.get_u() = x + *v.get_u();
01294     *z.get_u_x() = *v.get_u_x();
01295     *z.get_u_y() = *v.get_u_y();
01296     *z.get_u_z() = *v.get_u_z();
01297     *z.get_u_xx() = *v.get_u_xx();
01298     *z.get_u_xy() = *v.get_u_xy();
01299     *z.get_u_xz() = *v.get_u_xz();
01300     *z.get_u_yy() = *v.get_u_yy();
01301     *z.get_u_yz() = *v.get_u_yz();
01302     *z.get_u_zz() = *v.get_u_zz();
01303     *z.get_u_xxx() = *v.get_u_xxx();
01304     *z.get_u_xxy() = *v.get_u_xxy();
01305     *z.get_u_xxz() = *v.get_u_xxz();
01306     *z.get_u_xyy() = *v.get_u_xyy();
01307     *z.get_u_xyz() = *v.get_u_xyz();
01308     *z.get_u_xzz() = *v.get_u_xzz();
01309     *z.get_u_yyy() = *v.get_u_yyy();
01310     *z.get_u_yyz() = *v.get_u_yyz();
01311     *z.get_u_yzz() = *v.get_u_yzz();
01312     *z.get_u_zzz() = *v.get_u_zzz();
01313 
01314 
01315 
01316     return z;
01317   }
01318 
01323   df3_three_variable operator + (const df3_three_variable& x,
01324     const df3_three_variable& y)
01325   {
01326     df3_three_variable z;
01327     *z.get_u() = *x.get_u() + *y.get_u();
01328     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
01329     *z.get_u_z() = *x.get_u_z() + *y.get_u_z();
01330     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
01331     *z.get_u_xx() = *x.get_u_xx()+ *y.get_u_xx();
01332     *z.get_u_xy() = *x.get_u_xy()+ *y.get_u_xy();
01333     *z.get_u_xz() = *x.get_u_xz()+ *y.get_u_xz();
01334     *z.get_u_yy() = *x.get_u_yy()+ *y.get_u_yy();
01335     *z.get_u_yz() = *x.get_u_yz()+ *y.get_u_yz();
01336     *z.get_u_zz() = *x.get_u_zz()+ *y.get_u_zz();
01337     *z.get_u_xxx() = *x.get_u_xxx()+ *y.get_u_xxx();
01338     *z.get_u_xxy() = *x.get_u_xxy()+ *y.get_u_xxy();
01339     *z.get_u_xxz() = *x.get_u_xxz()+ *y.get_u_xxz();
01340     *z.get_u_xyy() = *x.get_u_xyy()+ *y.get_u_xyy();
01341     *z.get_u_xyz() = *x.get_u_xyz()+ *y.get_u_xyz();
01342     *z.get_u_xzz() = *x.get_u_xzz()+ *y.get_u_xzz();
01343     *z.get_u_yyy() = *x.get_u_yyy()+ *y.get_u_yyy();
01344     *z.get_u_yyz() = *x.get_u_yyz()+ *y.get_u_yyz();
01345     *z.get_u_yzz() = *x.get_u_yzz()+ *y.get_u_yzz();
01346     *z.get_u_zzz() = *x.get_u_zzz()+ *y.get_u_zzz();
01347     return z;
01348   }
01349 
01354   df3_three_variable operator - (const df3_three_variable& x,
01355     const df3_three_variable& y)
01356   {
01357     df3_three_variable z;
01358     *z.get_u() = *x.get_u() - *y.get_u();
01359     *z.get_u_x() = *x.get_u_x() - *y.get_u_x();
01360     *z.get_u_z() = *x.get_u_z() - *y.get_u_z();
01361     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
01362     *z.get_u_xx() = *x.get_u_xx() - *y.get_u_xx();
01363     *z.get_u_xy() = *x.get_u_xy() - *y.get_u_xy();
01364     *z.get_u_xz() = *x.get_u_xz() - *y.get_u_xz();
01365     *z.get_u_yy() = *x.get_u_yy() - *y.get_u_yy();
01366     *z.get_u_yz() = *x.get_u_yz() - *y.get_u_yz();
01367     *z.get_u_zz() = *x.get_u_zz() - *y.get_u_zz();
01368     *z.get_u_xxx() = *x.get_u_xxx() - *y.get_u_xxx();
01369     *z.get_u_xxy() = *x.get_u_xxy() - *y.get_u_xxy();
01370     *z.get_u_xxz() = *x.get_u_xxz() - *y.get_u_xxz();
01371     *z.get_u_xyy() = *x.get_u_xyy() - *y.get_u_xyy();
01372     *z.get_u_xyz() = *x.get_u_xyz() - *y.get_u_xyz();
01373     *z.get_u_xzz() = *x.get_u_xzz() - *y.get_u_xzz();
01374     *z.get_u_yyy() = *x.get_u_yyy() - *y.get_u_yyy();
01375     *z.get_u_yyz() = *x.get_u_yyz() - *y.get_u_yyz();
01376     *z.get_u_yzz() = *x.get_u_yzz() - *y.get_u_yzz();
01377     *z.get_u_zzz() = *x.get_u_zzz() - *y.get_u_zzz();
01378     return z;
01379   }
01380 
01385   df3_three_variable operator - (double x,
01386     const df3_three_variable& v)
01387   {
01388     df3_three_variable z;
01389 
01390     *z.get_u() = x - *v.get_u();
01391     *z.get_u_x() = -*v.get_u_x();
01392     *z.get_u_y() = -*v.get_u_y();
01393     *z.get_u_z() = -*v.get_u_z();
01394     *z.get_u_xx() = -*v.get_u_xx();
01395     *z.get_u_xy() = -*v.get_u_xy();
01396     *z.get_u_xz() = -*v.get_u_xz();
01397     *z.get_u_yy() = -*v.get_u_yy();
01398     *z.get_u_yz() = -*v.get_u_yz();
01399     *z.get_u_zz() = -*v.get_u_zz();
01400     *z.get_u_xxx() = -*v.get_u_xxx();
01401     *z.get_u_xxy() = -*v.get_u_xxy();
01402     *z.get_u_xxz() = -*v.get_u_xxz();
01403     *z.get_u_xyy() = -*v.get_u_xyy();
01404     *z.get_u_xyz() = -*v.get_u_xyz();
01405     *z.get_u_xzz() = -*v.get_u_xzz();
01406     *z.get_u_yyy() = -*v.get_u_yyy();
01407     *z.get_u_yyz() = -*v.get_u_yyz();
01408     *z.get_u_yzz() = -*v.get_u_yzz();
01409     *z.get_u_zzz() = -*v.get_u_zzz();
01410 
01411     return z;
01412   }
01413 
01418   df3_three_variable operator - (const df3_three_variable& v,
01419     double y)
01420   {
01421     df3_three_variable z;
01422 
01423     *z.get_u() =  *v.get_u()-y;
01424     *z.get_u_x() = *v.get_u_x();
01425     *z.get_u_y() = *v.get_u_y();
01426     *z.get_u_z() = *v.get_u_z();
01427     *z.get_u_xx() = *v.get_u_xx();
01428     *z.get_u_xy() = *v.get_u_xy();
01429     *z.get_u_xz() = *v.get_u_xz();
01430     *z.get_u_yy() = *v.get_u_yy();
01431     *z.get_u_yz() = *v.get_u_yz();
01432     *z.get_u_zz() = *v.get_u_zz();
01433     *z.get_u_xxx() = *v.get_u_xxx();
01434     *z.get_u_xxy() = *v.get_u_xxy();
01435     *z.get_u_xxz() = *v.get_u_xxz();
01436     *z.get_u_xyy() = *v.get_u_xyy();
01437     *z.get_u_xyz() = *v.get_u_xyz();
01438     *z.get_u_xzz() = *v.get_u_xzz();
01439     *z.get_u_yyy() = *v.get_u_yyy();
01440     *z.get_u_yyz() = *v.get_u_yyz();
01441     *z.get_u_yzz() = *v.get_u_yzz();
01442     *z.get_u_zzz() = *v.get_u_zzz();
01443 
01444     return z;
01445   }
01450   init_df3_three_variable::~init_df3_three_variable()
01451   {
01452     num_local_ind_var--;
01453     if (num_local_ind_var<0)
01454     {
01455       cerr << "This can't happen" << endl;
01456       ad_exit(1);
01457     }
01458     if (num_local_ind_var==0)
01459     {
01460       num_ind_var=0;
01461     }
01462   }
01463 
01468   init_df3_three_variable::init_df3_three_variable(const df1b2variable& _v)
01469   {
01470     ADUNCONST(df1b2variable,v)
01471     if (num_local_ind_var>2)
01472     {
01473       cerr << "can only have 3 independent_variables in df3_three_variable"
01474        " function" << endl;
01475       ad_exit(1);
01476     }
01477     ind_var[num_ind_var++]=&v;
01478     num_local_ind_var++;
01479     *get_u() =  *v.get_u();
01480     *get_u_x() = 0.0;
01481     *get_u_y() = 0.0;
01482     *get_u_z() = 0.0;
01483     switch(num_local_ind_var)
01484     {
01485     case 1:
01486       *get_u_x() = 1.0;
01487       break;
01488     case 2:
01489       *get_u_y() = 1.0;
01490       break;
01491     case 3:
01492       *get_u_z() = 1.0;
01493       break;
01494     default:
01495       cerr << "illegal num_ind_var value of " << num_ind_var
01496            << " in  df3_three_variable function" << endl;
01497       ad_exit(1);
01498     }
01499     *get_u_xx() = 0.0;
01500     *get_u_xy() = 0.0;
01501     *get_u_xz() = 0.0;
01502     *get_u_yy() = 0.0;
01503     *get_u_yz() = 0.0;
01504     *get_u_zz() = 0.0;
01505     *get_u_xxx() = 0.0;
01506     *get_u_xxy() = 0.0;
01507     *get_u_xxz() = 0.0;
01508     *get_u_xyy() = 0.0;
01509     *get_u_xyz() = 0.0;
01510     *get_u_xzz() = 0.0;
01511     *get_u_yyy() = 0.0;
01512     *get_u_yyz() = 0.0;
01513     *get_u_yzz() = 0.0;
01514     *get_u_zzz() = 0.0;
01515   }
01516 
01521   init_df3_three_variable::init_df3_three_variable(double v)
01522   {
01523     if (num_local_ind_var>2)
01524     {
01525       cerr << "can only have 3 independent_variables in df3_three_variable"
01526        " function" << endl;
01527       ad_exit(1);
01528     }
01529     num_local_ind_var++;
01530     *get_u() =  v;
01531     *get_u_x() = 0.0;
01532     *get_u_y() = 0.0;
01533     *get_u_z() = 0.0;
01534     switch(num_local_ind_var)
01535     {
01536     case 1:
01537       *get_u_x() = 1.0;
01538       break;
01539     case 2:
01540       *get_u_y() = 1.0;
01541       break;
01542     case 3:
01543       *get_u_z() = 1.0;
01544       break;
01545     default:
01546       cerr << "illegal num_ind_var value of " << num_ind_var
01547            << " in  df3_three_variable function" << endl;
01548       ad_exit(1);
01549     }
01550     *get_u_xx() = 0.0;
01551     *get_u_xy() = 0.0;
01552     *get_u_xz() = 0.0;
01553     *get_u_yy() = 0.0;
01554     *get_u_yz() = 0.0;
01555     *get_u_zz() = 0.0;
01556     *get_u_xxx() = 0.0;
01557     *get_u_xxy() = 0.0;
01558     *get_u_xxz() = 0.0;
01559     *get_u_xyy() = 0.0;
01560     *get_u_xyz() = 0.0;
01561     *get_u_xzz() = 0.0;
01562     *get_u_yyy() = 0.0;
01563     *get_u_yyz() = 0.0;
01564     *get_u_yzz() = 0.0;
01565     *get_u_zzz() = 0.0;
01566   }
01567 
01568   df3_three_variable::df3_three_variable(void)
01569   {
01570   }
01571 
01576 df3_three_matrix choleski_decomp(const df3_three_matrix& MM)
01577 {
01578   // kludge to deal with constantness
01579   df3_three_matrix & M= (df3_three_matrix &) MM;
01580   int rmin=M.indexmin();
01581   int cmin=M(rmin).indexmin();
01582   int rmax=M.indexmax();
01583   int cmax=M(rmin).indexmax();
01584   if (rmin !=1 || cmin !=1)
01585   {
01586     cerr << "minimum row and column inidices must equal 1 in "
01587       "df1b2matrix choleski_decomp(const df3_three_atrix& MM)"
01588          << endl;
01589     ad_exit(1);
01590   }
01591   if (rmax !=cmax)
01592   {
01593     cerr << "Error in df1b2matrix choleski_decomp(const df3_three_matrix& MM)"
01594       " Matrix not square" << endl;
01595     ad_exit(1);
01596   }
01597 
01598   int n=rmax-rmin+1;
01599   df3_three_matrix L(1,n,1,n);
01600 #ifndef SAFE_INITIALIZE
01601     L.initialize();
01602 #endif
01603 
01604   int i,j,k;
01605   df3_three_variable tmp;
01606 
01607     if (value(M(1,1))<=0)
01608     {
01609       cerr << "Error matrix not positive definite in choleski_decomp"
01610         <<endl;
01611       ad_exit(1);
01612     }
01613 
01614   L(1,1)=sqrt(M(1,1));
01615   for (i=2;i<=n;i++)
01616   {
01617     L(i,1)=M(i,1)/L(1,1);
01618   }
01619 
01620   for (i=2;i<=n;i++)
01621   {
01622     for (j=2;j<=i-1;j++)
01623     {
01624       tmp=M(i,j);
01625       for (k=1;k<=j-1;k++)
01626       {
01627         tmp-=L(i,k)*L(j,k);
01628       }
01629       L(i,j)=tmp/L(j,j);
01630     }
01631     tmp=M(i,i);
01632     for (k=1;k<=i-1;k++)
01633     {
01634       tmp-=L(i,k)*L(i,k);
01635     }
01636 
01637     if (value(tmp)<=0)
01638     {
01639       cerr << "Error matrix not positive definite in choleski_decomp"
01640         <<endl;
01641       ad_exit(1);
01642     }
01643 
01644     L(i,i)=sqrt(tmp);
01645   }
01646 
01647   return L;
01648 }