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