ADMB Documentation  11.1.2503
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f12.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f12.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00013 void ad_read_pass2(void);
00014 
00019  int df1b2_gradlist::write_pass1(const df1b2variable * _px,
00020    const df1b2variable * _py,df1b2variable * pz,df1b2function2 * pf)
00021  {
00022    ADUNCONST(df1b2variable*,px)
00023    ADUNCONST(df1b2variable*,py)
00024    ncount++;
00025 #if defined(CHECK_COUNT)
00026   if (ncount >= ncount_check)
00027     ncount_checker(ncount,ncount_check);
00028 #endif
00029    int nvar=df1b2variable::nvar;
00030 
00031    int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00032      +2*(nvar+1)*sizeof(double);
00033 // string identifier debug stuff
00034 #if defined(SAFE_ALL)
00035   char ids[]="BY";
00036   int slen=strlen(ids);
00037   total_bytes+=slen;
00038 #endif
00039   list.check_buffer_size(total_bytes);
00040   void * tmpptr=list.bptr;
00041 #if defined(SAFE_ALL)
00042   memcpy(list,ids,slen);
00043 #endif
00044 // end of string identifier debug stuff
00045 
00046    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00047    memcpy(list,(df1b2_header*)(py),sizeof(df1b2_header));
00048    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00049    memcpy(list,&pf,sizeof(char *));
00050    //*(char**)(list.bptr)=(char*)pf;
00051    memcpy(list,px->get_u(),sizeof(double));
00052    memcpy(list,py->get_u(),sizeof(double));
00053    memcpy(list,px->get_u_dot(),nvar*sizeof(double));
00054    memcpy(list,py->get_u_dot(),nvar*sizeof(double));
00055    // ***** write  record size
00056    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00057    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2);
00058    ++nlist;
00059    return 0;
00060  }
00061 
00062 
00063 void read_pass2_1(void);
00064 void read_pass2_2(void);
00065 void read_pass2_3(void);
00066 
00071 void ad_read_pass2(void)
00072 {
00073   switch(df1b2variable::passnumber)
00074   {
00075   case 1:
00076     read_pass2_1();
00077     break;
00078   case 2:
00079     read_pass2_2();
00080     break;
00081   case 3:
00082     read_pass2_3();
00083     break;
00084   default:
00085     cerr << "illegal value for df1b2variable::pass = "
00086          << df1b2variable::passnumber << endl;
00087     exit(1);
00088   }
00089 }
00090 
00095 void read_pass2_1(void)
00096 {
00097   // We are going backword for bptr and nbptr
00098   // and  forward for bptr2 and nbptr2
00099   // the current entry+2 in bptr is the size of the record i.e
00100   // points to the next record
00101   //char * bptr=f1b2gradlist->bptr;
00102   //char * bptr2=f1b2gradlist2->bptr;
00103   int nvar=df1b2variable::nvar;
00104   test_smartlist& list=f1b2gradlist->list;
00105   //f1b2gradlist->nlist-=sizeof(int);
00106   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00107   list-=num_bytes;
00108   list.saveposition(); // save pointer to beginning of record;
00109   double xu,yu;
00110   //ad_dstar xdot,ydot;
00111   df1b2function2 * pf;
00112 
00113   // get info from tape1
00114 #if defined(SAFE_ARRAYS)
00115   checkidentiferstring("BY",f1b2gradlist->list);
00116 #endif
00117   char * bptr=f1b2gradlist->list.bptr;
00118   df1b2_header * px=(df1b2_header *) bptr;
00119   bptr+=sizeof(df1b2_header);
00120   df1b2_header * py=(df1b2_header *) bptr;
00121   bptr+=sizeof(df1b2_header);
00122   df1b2_header * pz=(df1b2_header *) bptr;
00123   bptr+=sizeof(df1b2_header);
00124   pf=*(df1b2function2 **) bptr;
00125   bptr+=sizeof(char*);
00126   memcpy(&xu,bptr,sizeof(double));
00127   bptr+=sizeof(double);
00128   memcpy(&yu,bptr,sizeof(double));
00129   bptr+=sizeof(double);
00130   double * xdot=(double*)bptr;
00131   bptr+=nvar*sizeof(double);
00132   double * ydot=(double*)bptr;
00133 
00134   list.restoreposition(); // save pointer to beginning of record;
00135 
00136   // ****************************************************************
00137   // turn this off if no third derivatives are calculated
00138   // if (!no_third_derivatives)
00139   // {
00140   // save for second reverse pass
00141   // save identifier 1
00142      test_smartlist & list2 = f1b2gradlist->list2;
00143 
00144 
00145    int total_bytes=2*nvar*sizeof(double);
00146 // string identifier debug stuff
00147 #if defined(SAFE_ALL)
00148   char ids[]="FW";
00149   int slen=strlen(ids);
00150   total_bytes+=slen;
00151 #endif
00152   list2.check_buffer_size(total_bytes);
00153   void * tmpptr=list2.bptr;
00154 #if defined(SAFE_ALL)
00155   memcpy(list2,ids,slen);
00156 #endif
00157 
00158      fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00159      memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00160      memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00161      *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00162      ++nlist2;
00163   // }
00164   //
00165   // ****************************************************************
00166 #if defined(PRINT_DERS)
00167  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00168   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00169   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00170   (pf->d3f222)(xu,yu),1);
00171  print_derivatives(pz,"z");
00172  print_derivatives(px,"x");
00173  print_derivatives(py,"y");
00174 #endif
00175 #if defined(__DERCHECK__)
00176   if (derchecker)
00177   if (derchecker->node_number)
00178   {
00179     if (derchecker->counter == derchecker->node_number)
00180     {
00181       switch (derchecker->pass_number) // increment the variable of interest
00182       {
00183       case 2:
00184         switch(derchecker->vartype)
00185         {
00186         case 1:
00187           if (!derchecker->dotflag)
00188             px->u_bar[derchecker->index-1]+=derchecker->delta;
00189           else
00190             px->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00191           break;
00192         case 2:
00193           if (!derchecker->dotflag)
00194             py->u_bar[derchecker->index-1]+=derchecker->delta;
00195           else
00196             py->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00197           break;
00198         case 3:
00199           if (!derchecker->dotflag)
00200             pz->u_bar[derchecker->index-1]+=derchecker->delta;
00201           else
00202             pz->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00203           break;
00204         default:
00205           cerr << "Invalid index value for dercheck_index was "
00206                << derchecker->index << endl;
00207           break;
00208         }
00209         break;
00210       case 3:
00211         switch(derchecker->vartype)
00212         {
00213         case 1:
00214           if (!derchecker->dotflag)
00215             px->u_bar[derchecker->index-1]-=derchecker->delta;
00216           else
00217             px->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00218           break;
00219         case 2:
00220           if (!derchecker->dotflag)
00221             py->u_bar[derchecker->index-1]-=derchecker->delta;
00222           else
00223             py->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00224           break;
00225         case 3:
00226           if (!derchecker->dotflag)
00227             pz->u_bar[derchecker->index-1]-=derchecker->delta;
00228           else
00229             pz->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00230           break;
00231         default:
00232           cerr << "Invalid index value for dercheck_index was "
00233                << derchecker->index << endl;
00234           break;
00235         }
00236         break;
00237       }
00238     }
00239   }
00240 #endif
00241 
00242   // Do first reverse pass calculations
00243   int i;
00244   for (i=0;i<nvar;i++)
00245   {
00246     //double df1=(pf->df1)(xu,yu);
00247     px->u_bar[i]+=(pf->df1)(xu,yu)*pz->u_bar[i];
00248   }
00249   for (i=0;i<nvar;i++)
00250   {
00251     py->u_bar[i]+=(pf->df2)(xu,yu)*pz->u_bar[i];
00252   }
00253   for (i=0;i<nvar;i++)
00254   {
00255     px->u_bar[i]+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00256     px->u_bar[i]+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00257 #if defined(ADDEBUG_PRINT)
00258     cout << px->u_bar[i] << " " << pz->u_dot_bar[i] << " " << addebug_count
00259          << endl;
00260 #endif
00261   }
00262   for (i=0;i<nvar;i++)
00263   {
00264     //py->u_bar[i]+=(pf->d2f22)(*(px->u),*(py->u))*ydot[i]*pz->u_dot_bar[i];
00265     //py->u_bar[i]+=(pf->d2f12)(*(px->u),*(py->u))*xdot[i]*pz->u_dot_bar[i];
00266     py->u_bar[i]+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_bar[i];
00267     py->u_bar[i]+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00268 #if defined(ADDEBUG_PRINT)
00269     cout << py->u_bar[i] << " " << pz->u_dot_bar[i] << " " << addebug_count
00270          << endl;
00271 #endif
00272   }
00273   for (i=0;i<nvar;i++)
00274   {
00275     //px->u_dot_bar[i]+=(pf->df1)(*(px->u),*(py->u))*pz->u_dot_bar[i];
00276     px->u_dot_bar[i]+=(pf->df1)(xu,yu)*pz->u_dot_bar[i];
00277 #if defined(ADDEBUG_PRINT)
00278     cout << px->u_dot_bar[i] << " " << addebug_count << endl;
00279     cout << pz->u_dot_bar[i] << " " << addebug_count << endl;
00280 #endif
00281   }
00282   for (i=0;i<nvar;i++)
00283   {
00284     py->u_dot_bar[i]+=(pf->df2)(xu,yu)*pz->u_dot_bar[i];
00285 #if defined(ADDEBUG_PRINT)
00286     cout << py->u_dot_bar[i] << " " << addebug_count << endl;
00287     cout << pz->u_dot_bar[i] << " " << addebug_count << endl;
00288 #endif
00289   }
00290 
00291   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00292   for (i=0;i<nvar;i++)
00293   {
00294     pz->u_bar[i]=0;
00295   }
00296   for (i=0;i<nvar;i++)
00297   {
00298     pz->u_dot_bar[i]=0;
00299   }
00300 
00301 #if defined(PRINT_DERS)
00302  print_derivatives(pz,"z");
00303  print_derivatives(px,"x");
00304  print_derivatives(py,"y");
00305 #endif
00306 }
00307 
00312 void read_pass2_2(void)
00313 {
00314   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00315   // We are going forward for bptr and backword for bptr2
00316   //
00317   // list 1
00318   //
00319   int nvar=df1b2variable::nvar;
00320   test_smartlist & list=f1b2gradlist->list;
00321 
00322   int total_bytes=3*sizeof(df1b2_header)+sizeof(char*)
00323     +2*(nvar+1)*sizeof(double);
00324 // string identifier debug stuff
00325 #if defined(SAFE_ALL)
00326   char ids[]="BY";
00327   int slen=strlen(ids);
00328   total_bytes+=slen;
00329 #endif
00330   list.check_buffer_size(total_bytes);
00331 // end of string identifier debug stuff
00332 
00333   list.saveposition(); // save pointer to beginning of record;
00334   fixed_smartlist & nlist=f1b2gradlist->nlist;
00335    // nlist-=sizeof(int);
00336   // get record size
00337   int num_bytes=nlist.bptr->numbytes;
00338     // nlist+=nlist_record_size;
00339   //
00340   // list 2
00341   //
00342   test_smartlist & list2=f1b2gradlist->list2;
00343   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00344   // get record size
00345   int num_bytes2=*nlist2.bptr;
00346   --nlist2;
00347   // backup the size of the record
00348   list2-=num_bytes2;
00349   list2.saveposition(); // save pointer to beginning of record;
00350   // save the pointer to the beginning of the record
00351   // bptr and bptr2 now both point to the beginning of their records
00352 
00353   double xu,yu;
00354   double * xdot;
00355   double * ydot;
00356   //df1b2_header x,z;
00357   df1b2function2 * pf;
00358 
00359   // get info from tape1
00360   // get info from tape1
00361 #if defined(SAFE_ARRAYS)
00362   checkidentiferstring("BY",list);
00363   checkidentiferstring("FW",list2);
00364 #endif
00365   df1b2_header * px=(df1b2_header *) list.bptr;
00366   list.bptr+=sizeof(df1b2_header);
00367   df1b2_header * py=(df1b2_header *) list.bptr;
00368   list.bptr+=sizeof(df1b2_header);
00369   df1b2_header * pz=(df1b2_header *) list.bptr;
00370   list.bptr+=sizeof(df1b2_header);
00371   pf=*(df1b2function2 **) list.bptr;
00372   list.bptr+=sizeof(char*);
00373   memcpy(&xu,list.bptr,sizeof(double));
00374   list.bptr+=sizeof(double);
00375   memcpy(&yu,list.bptr,sizeof(double));
00376   list.bptr+=sizeof(double);
00377   xdot=(double*)list.bptr;
00378   list.bptr+=nvar*sizeof(double);
00379   ydot=(double*)list.bptr;
00380   list.restoreposition(num_bytes); // save pointer to beginning of record;
00381 
00382   double * zbar;
00383   double * zdotbar;
00384 
00385   zbar=(double*)list2.bptr;
00386   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00387   list2.restoreposition(); // save pointer to beginning of record;
00388 
00389   double * x_tilde=px->get_u_tilde();
00390   double * x_dot_tilde=px->get_u_dot_tilde();
00391   double * x_bar_tilde=px->get_u_bar_tilde();
00392   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00393   double * y_tilde=py->get_u_tilde();
00394   double * y_dot_tilde=py->get_u_dot_tilde();
00395   double * y_bar_tilde=py->get_u_bar_tilde();
00396   double * y_dot_bar_tilde=py->get_u_dot_bar_tilde();
00397   double * z_bar_tilde=pz->get_u_bar_tilde();
00398   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00399   // Do second "reverse-reverse" pass calculations
00400 #if defined(PRINT_DERS)
00401  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00402   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00403   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00404   (pf->d3f222)(xu,yu),1);
00405  print_derivatives(pz,"z");
00406  print_derivatives(px,"x");
00407  print_derivatives(py,"y");
00408 #endif
00409 
00410   int i;
00411 
00412   for (i=0;i<nvar;i++)
00413   {
00414     z_bar_tilde[i]=0;
00415     z_dot_bar_tilde[i]=0;
00416   }
00417 
00418   // start with x and add y
00419   for (i=0;i<nvar;i++)
00420   {
00421     *x_tilde+=(pf->d2f11)(xu,yu)*zbar[i]*x_bar_tilde[i];
00422     z_bar_tilde[i]+=(pf->df1)(xu,yu)*x_bar_tilde[i];
00423     *y_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*x_bar_tilde[i];
00424   }
00425 
00426   for (i=0;i<nvar;i++)
00427   {
00428     *x_tilde+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00429     *y_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00430     z_dot_bar_tilde[i]+=(pf->df1)(xu,yu)*x_dot_bar_tilde[i];
00431   }
00432 
00433   for (i=0;i<nvar;i++)
00434   {
00435     x_dot_tilde[i]+=(pf->d2f11)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00436     z_dot_bar_tilde[i]+=(pf->d2f11)(xu,yu)*xdot[i]*x_bar_tilde[i];
00437     *x_tilde+=(pf->d3f111)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00438     *y_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00439   }
00440   // start with y and add x
00441   for (i=0;i<nvar;i++)
00442   {
00443     *y_tilde+=(pf->d2f22)(xu,yu)*zbar[i]*y_bar_tilde[i];
00444     *x_tilde+=(pf->d2f12)(xu,yu)*zbar[i]*y_bar_tilde[i];
00445     z_bar_tilde[i]+=(pf->df2)(xu,yu)*y_bar_tilde[i];
00446   }
00447 
00448   for (i=0;i<nvar;i++)
00449   {
00450     *y_tilde+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00451     *x_tilde+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_dot_bar_tilde[i];
00452     z_dot_bar_tilde[i]+=(pf->df2)(xu,yu)*y_dot_bar_tilde[i];
00453   }
00454 
00455   for (i=0;i<nvar;i++)
00456   {
00457     y_dot_tilde[i]+=(pf->d2f22)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00458     z_dot_bar_tilde[i]+=(pf->d2f22)(xu,yu)*ydot[i]*y_bar_tilde[i];
00459     *y_tilde+=(pf->d3f222)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00460     *x_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*y_bar_tilde[i];
00461   }
00462   for (i=0;i<nvar;i++)
00463   {
00464     *x_tilde+=(pf->d3f112)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00465     *y_tilde+=(pf->d3f122)(xu,yu)*ydot[i]*zdotbar[i]*x_bar_tilde[i];
00466     y_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00467     z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*ydot[i]*x_bar_tilde[i];
00468   }
00469   for (i=0;i<nvar;i++)
00470   {
00471     *x_tilde+=(pf->d3f112)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00472     *y_tilde+=(pf->d3f122)(xu,yu)*xdot[i]*zdotbar[i]*y_bar_tilde[i];
00473     x_dot_tilde[i]+=(pf->d2f12)(xu,yu)*zdotbar[i]*y_bar_tilde[i];
00474     z_dot_bar_tilde[i]+=(pf->d2f12)(xu,yu)*xdot[i]*y_bar_tilde[i];
00475   }
00476 #if defined(__DERCHECK__)
00477   if (derchecker->node_number)
00478   {
00479     if (derchecker->counter == derchecker->node_number)
00480     {
00481       if (derchecker->pass_number==1) // increment the variable of interest
00482       {
00483         switch(derchecker->vartype)
00484         {
00485         case 1:
00486           if (!derchecker->dotflag)
00487             derchecker->der_value=
00488               px->u_bar_tilde[derchecker->index-1];
00489           else
00490             derchecker->der_value=
00491               px->u_dot_bar_tilde[derchecker->index-1];
00492           break;
00493         case 2:
00494           if (!derchecker->dotflag)
00495             derchecker->der_value=
00496               py->u_bar_tilde[derchecker->index-1];
00497           else
00498             derchecker->der_value=
00499               py->u_dot_bar_tilde[derchecker->index-1];
00500           break;
00501         case 3:
00502           if (!derchecker->dotflag)
00503             derchecker->der_value=
00504               pz->u_bar_tilde[derchecker->index-1];
00505           else
00506             derchecker->der_value=
00507               pz->u_dot_bar_tilde[derchecker->index-1];
00508           break;
00509         default:
00510           cerr << "Invalid index value for dercheck_index was "
00511                << derchecker->index << endl;
00512         }
00513       }
00514     }
00515   }
00516 #endif
00517 #if defined(PRINT_DERS)
00518  print_derivatives(pz,"z");
00519  print_derivatives(px,"x");
00520  print_derivatives(py,"y");
00521 #endif
00522 }
00523 
00528 void read_pass2_3(void)
00529 {
00530   // We are going backword for bptr and forward for bptr2
00531   // the current entry+2 in bptr is the size of the record i.e
00532   // points to the next record
00533   int nvar=df1b2variable::nvar;
00534   fixed_smartlist & nlist=f1b2gradlist->nlist;
00535   test_smartlist& list=f1b2gradlist->list;
00536    // nlist-=sizeof(int);
00537   // get record size
00538   int num_bytes=nlist.bptr->numbytes;
00539   // backup the size of the record
00540   list-=num_bytes;
00541   list.saveposition(); // save pointer to beginning of record;
00542   // save the pointer to the beginning of the record
00543   double xu;
00544   double yu;
00545   double * xdot;
00546   double * ydot;
00547   //df1b2_header x,z;
00548   df1b2function2 * pf;
00549 
00550   // get info from tape1
00551   // get info from tape1
00552 #if defined(SAFE_ARRAYS)
00553   checkidentiferstring("BY",list);
00554 #endif
00555   df1b2_header * px=(df1b2_header *) list.bptr;
00556   list.bptr+=sizeof(df1b2_header);
00557   df1b2_header * py=(df1b2_header *) list.bptr;
00558   list.bptr+=sizeof(df1b2_header);
00559   df1b2_header * pz=(df1b2_header *) list.bptr;
00560   list.bptr+=sizeof(df1b2_header);
00561   pf=*(df1b2function2 **) list.bptr;
00562   list.bptr+=sizeof(char*);
00563   memcpy(&xu,list.bptr,sizeof(double));
00564   list.bptr+=sizeof(double);
00565   memcpy(&yu,list.bptr,sizeof(double));
00566   list.bptr+=sizeof(double);
00567   xdot=(double*)list.bptr;
00568   list.bptr+=nvar*sizeof(double);
00569   ydot=(double*)list.bptr;
00570   list.restoreposition(); // save pointer to beginning of record;
00571   int i;
00572 #if defined(PRINT_DERS)
00573  print_derivatives(pf->funname,(pf->f)(xu,yu),(pf->df1)(xu,yu),
00574   (pf->df2)(xu,yu),(pf->d2f11)(xu,yu),(pf->d2f12)(xu,yu),(pf->d2f22)(xu,yu),
00575   (pf->d3f111)(xu,yu),(pf->d3f112)(xu,yu),(pf->d3f122)(xu,yu),
00576   (pf->d3f222)(xu,yu),1);
00577  print_derivatives(pz,"z");
00578  print_derivatives(px,"x");
00579  print_derivatives(py,"y");
00580 #endif
00581 
00582   *(px->u_tilde)+=(pf->df1)(xu,yu)* *(pz->u_tilde);
00583   *(py->u_tilde)+=(pf->df2)(xu,yu)* *(pz->u_tilde);
00584   for (i=0;i<nvar;i++)
00585   {
00586     *(px->u_tilde)+=(pf->d2f11)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00587     *(py->u_tilde)+=(pf->d2f12)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00588     *(py->u_tilde)+=(pf->d2f22)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00589     *(px->u_tilde)+=(pf->d2f12)(xu,yu)*ydot[i]*pz->u_dot_tilde[i];
00590   }
00591   for (i=0;i<nvar;i++)
00592   {
00593     px->u_dot_tilde[i]+=(pf->df1)(xu,yu)*pz->u_dot_tilde[i];
00594     py->u_dot_tilde[i]+=(pf->df2)(xu,yu)*pz->u_dot_tilde[i];
00595   }
00596   *(pz->u_tilde)=0;
00597   for (i=0;i<nvar;i++)
00598   {
00599     pz->u_dot_tilde[i]=0;
00600   }
00601 #if defined(PRINT_DERS)
00602  print_derivatives(pz,"z");
00603  print_derivatives(px,"x");
00604  print_derivatives(py,"y");
00605 #endif
00606 }