source: trunk/ADOL-C/examples/additional_examples/timing/vfunc_shuttle.cpp @ 42

Last change on this file since 42 was 42, checked in by awalther, 10 years ago

set svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 7.6 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     vfunc_shuttle.cpp
4 Revision: $Id: vfunc_shuttle.cpp 42 2009-07-15 18:37:17Z awalther $
5 Contents: Example of function module containing the shuttle example
6            (based on shuttlexam.c of version 1.7)
7 
8   Each << function module >> contains:
9         
10     (1) const char* const controlFileName
11     (2) int indepDim;
12     (3) int depDim;
13     (4) void initProblemParameters( void )
14     (5) void initIndependents( double* indeps )
15     (6) void originalVectorFunction( double* indeps, double* deps )
16     (7) void tapingVectorFunction( int tag, double* indeps, double* deps )   
17 
18 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
19               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
20   
21 This file is part of ADOL-C. This software is provided as open source.
22 Any use, reproduction, or distribution of the software constitutes
23 recipient's acceptance of the terms of the accompanying license file.
24 
25---------------------------------------------------------------------------*/
26#define _VFUNC_SHUTTLE_C_
27
28
29/****************************************************************************/
30/*                                                                 INCLUDES */
31#include <adolc.h>
32
33#include <cmath>
34
35
36/****************************************************************************/
37/*                                                         GLOBAL VARIABLES */
38
39/*--------------------------------------------------------------------------*/
40/*                                                        Control file name */
41const char* controlFileName = "shuttlexam.ctrl";
42
43/*--------------------------------------------------------------------------*/
44/*                                                               Dimensions */
45int indepDim;
46int depDim;
47
48/*--------------------------------------------------------------------------*/
49/*                                       Other problem dependent parameters */
50const double Pi = 3.141592654;
51const double ae = 20902900.0;
52const double mu = 0.14E+17;
53const double a  = 40.0;
54const double S  = 2690.0;
55const double crtd = 180.0/Pi;
56const double cl = 0.84-0.48*(38.0-a*crtd)/26.0;
57const double C0 =  3.974960446019;
58const double C1 = -0.01448947694635;
59const double C2 = -0.2156171551995e-4;
60const double C3 = -0.1089609507291e-7;
61const double V0 =  0.0;
62const double ma = 5964.496499824;
63const double Om = .72921159e-4;
64
65
66/****************************************************************************/
67/*                                                  INIT PROBLEM PARAMETERS */
68void initProblemParameters( void ) {
69    fprintf(stdout,"SHUTTLEXAM (ADOL-C Example)\n\n");
70
71    /* number of indeps & deps */
72    indepDim = 14;
73    depDim   = 7;
74}
75
76
77/****************************************************************************/
78/*                                                        INITIALIZE INDEPs */
79void initIndependents( double* indeps ) {
80    indeps[0]  = 264039.328;   /* H */
81    indeps[1]  = 177.718047;   /* x */
82    indeps[2]  = 32.0417885;   /* l */
83    indeps[3]  = 24317.0798;   /* V */
84    indeps[4]  = -0.749986488; /* g */
85    indeps[5]  = 62.7883367;   /* A */
86    indeps[6]  = 41.100771834; /* b */
87    indeps[7]  = -318;         /* Hp */
88    indeps[8]  = 0.01;         /* xp */
89    indeps[9]  = 0.1;          /* lp */
90    indeps[10] = -3.6;         /* Vp */
91    indeps[11] = 0.001;        /* gp */
92    indeps[12] = 0.1;          /* Ap */
93    indeps[13] = 0.06;         /* bp */
94}
95
96
97/****************************************************************************/
98/*                                                 ORIGINAL SCALAR FUNCTION */
99
100/*--------------------------------------------------------------------------*/
101/*                                                     The shuttle function */
102void shuttle( double* indeps, double* deps ) {
103    double r,gr,rho,L,cd,Z;
104    double sing,cosg,sinA,cosA,sinl,cosl,tanl;
105
106    r  = indeps[0]+ae;
107    gr = mu/(r*r);
108    rho= 0.002378*exp(-indeps[0]/23800.0);
109    L  = 0.5*rho*cl*S*indeps[3]*indeps[3];
110    cd = 0.78-0.58*(38.0-a*crtd)/26.0;
111    Z  = .5*rho*cd*S*indeps[3]*indeps[3];
112    // evaluate the dynamic equations ...
113    sinA = sin(indeps[5]);
114    cosA = cos(indeps[5]);
115    sing = sin(indeps[4]);
116    cosg = cos(indeps[4]);
117    sinl = sin(indeps[2]);
118    cosl = cos(indeps[2]);
119    tanl = sinl/cosl;
120    deps[0] = indeps[3]*sing-indeps[7];
121    deps[1] = indeps[3]*cosg*sinA/(r*cosl)-indeps[8];
122    deps[2] = indeps[3]*cosg*cosA/r-indeps[9];
123    deps[3] = -Z/ma-gr*sing-Om*Om*r*cosl
124              *(sinl*cosA*cosg-cosl*sing)-indeps[10];
125    deps[4] = L*cos(indeps[6])/(ma*indeps[3])
126              +cosl/indeps[3]*(indeps[3]*indeps[3]/r-gr)
127              +2*Om*cosl*sinA
128              +Om*Om*r*cosl/indeps[3]*(sinl*cosA*sing+cosl*cosg)
129              -indeps[11];
130    deps[5] = L*sin(indeps[6])/(ma*indeps[3]*cosg)+indeps[3]/r*cosg*sinA*tanl
131              -2*Om*(cosl*cosA*sing/cosg-sinl)
132              +Om*Om*r*cosl*sinl*sinA/(indeps[3]*cosg)-indeps[12];
133    deps[6] = Z/ma
134              -(C0+(indeps[3]-V0)*(C1+(indeps[3]-V0)*(C2+(indeps[3]-V0)*C3)));
135}
136
137/*--------------------------------------------------------------------------*/
138/*                                                   The interface function */
139void originalVectorFunction( double* indeps, double* deps ) {
140    shuttle(indeps,deps);
141}
142
143
144/****************************************************************************/
145/*                                                   TAPING SCALAR FUNCTION */
146
147/*--------------------------------------------------------------------------*/
148/*                                              The active shuttle function */
149void activeShuttle( adouble* indeps, adouble* deps ) {
150    adouble r,gr,rho,L,cd,Z;
151    adouble sing,cosg,sinA,cosA,sinl,cosl,tanl;
152
153    r  = indeps[0]+ae;
154    gr = mu/(r*r);
155    rho= 0.002378*exp(-indeps[0]/23800.0);
156    L  = 0.5*rho*cl*S*indeps[3]*indeps[3];
157    cd = 0.78-0.58*(38.0-a*crtd)/26.0;
158    Z  = .5*rho*cd*S*indeps[3]*indeps[3];
159    // evaluate the dynamic equations ...
160    sinA = sin(indeps[5]);
161    cosA = cos(indeps[5]);
162    sing = sin(indeps[4]);
163    cosg = cos(indeps[4]);
164    sinl = sin(indeps[2]);
165    cosl = cos(indeps[2]);
166    tanl = sinl/cosl;
167    deps[0] = indeps[3]*sing-indeps[7];
168    deps[1] = indeps[3]*cosg*sinA/(r*cosl)-indeps[8];
169    deps[2] = indeps[3]*cosg*cosA/r-indeps[9];
170    deps[3] = -Z/ma-gr*sing-Om*Om*r*cosl
171              *(sinl*cosA*cosg-cosl*sing)-indeps[10];
172    deps[4] = L*cos(indeps[6])/(ma*indeps[3])
173              +cosl/indeps[3]*(indeps[3]*indeps[3]/r-gr)
174              +2*Om*cosl*sinA
175              +Om*Om*r*cosl/indeps[3]*(sinl*cosA*sing+cosl*cosg)
176              -indeps[11];
177    deps[5] = L*sin(indeps[6])/(ma*indeps[3]*cosg)+indeps[3]/r*cosg*sinA*tanl
178              -2*Om*(cosl*cosA*sing/cosg-sinl)
179              +Om*Om*r*cosl*sinl*sinA/(indeps[3]*cosg)-indeps[12];
180    deps[6] = Z/ma
181              -(C0+(indeps[3]-V0)*(C1+(indeps[3]-V0)*(C2+(indeps[3]-V0)*C3)));
182}
183
184
185/*--------------------------------------------------------------------------*/
186/*                                                   The interface function */
187void tapingVectorFunction( int tag, double* indeps, double* deps ) {
188    int i;
189    trace_on(tag);
190    adouble* activeIndeps = new adouble[indepDim];
191    adouble* activeDeps   = new adouble[depDim];
192    adouble* aIP = activeIndeps;
193    double*  iP  = indeps;
194    for (i=0; i<indepDim; i++)
195        *aIP++ <<= *iP++;
196    activeShuttle(activeIndeps,activeDeps);
197    aIP = activeDeps;
198    iP  = deps;
199    for (i=0; i<depDim; i++)
200        *aIP++ >>= *iP++;
201    trace_off();
202}
203
204#undef _VFUNC_SHUTTLE_C_
205
206
207
208
209
Note: See TracBrowser for help on using the repository browser.