← Back to branch summary

~ubuntu-branches/ubuntu/vivid/gyoto/vivid

« back to all changes in this revision

Viewing changes to yorick/stdplug/gyoto_StarTrace.C

  • Committer: Package Import Robot
  • Author(s): Thibaut Paumard
  • Date: 2014-10-21 14:21:10 UTC
  • mfrom: (8.2.4 sid)
  • Revision ID: package-import@ubuntu.com-20141021142110-1u2odjj2r0hvdh5g
Tags: 0.2.3-1
* New upstream release
* Move to debian-astro team

Show diffs side-by-side

added added

removed removed

Lines of Context:
yorick/stdplug/gyoto_StarTrace.C
 
1
/*
 
2
    Copyright 2013 Thibaut Paumard
 
3
 
 
4
    This file is part of Gyoto.
 
5
 
 
6
    Gyoto is free software: you can redistribute it and/or modify
 
7
    it under the terms of the GNU General Public License as published by
 
8
    the Free Software Foundation, either version 3 of the License, or
 
9
    (at your option) any later version.
 
10
 
 
11
    Gyoto is distributed in the hope that it will be useful,
 
12
    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
    GNU General Public License for more details.
 
15
 
 
16
    You should have received a copy of the GNU General Public License
 
17
    along with Gyoto.  If not, see <http://www.gnu.org/licenses/>.
 
18
 */
 
19
 
 
20
#include "../ygyoto.h"
 
21
#include "yapi.h"
 
22
#include <GyotoStarTrace.h>
 
23
#ifdef GYOTO_USE_XERCES
 
24
#include <GyotoFactory.h>
 
25
#endif
 
26
 
 
27
#include <iostream>
 
28
#include <cstring>
 
29
using namespace std;
 
30
 
 
31
using namespace Gyoto;
 
32
using namespace Gyoto::Astrobj;
 
33
 
 
34
#define OBJ ao
 
35
 
 
36
// on_eval worker
 
37
void ygyoto_StarTrace_eval(SmartPointer<Astrobj::Generic>* ao_, int argc) {
 
38
  GYOTO_DEBUG << endl;
 
39
 
 
40
  // Define keywords
 
41
  static char const * knames[]={
 
42
    "unit",
 
43
    "radius", "metric", "initcoord", "spectrum", "opacity", "delta", "adaptive",
 
44
    "deltamaxoverradius", "deltamaxoverdistance", "tmin", "tmax",
 
45
    "maxiter", "reset", "xfill",
 
46
    YGYOTO_ASTROBJ_GENERIC_KW,
 
47
    "get_skypos", "get_txyz", "get_prime", "get_coord", "get_cartesian",
 
48
    "star",
 
49
    0
 
50
  };
 
51
 
 
52
  YGYOTO_WORKER_INIT(Astrobj, StarTrace, knames, YGYOTO_ASTROBJ_GENERIC_KW_N+21);
 
53
 
 
54
  YGYOTO_WORKER_SET_UNIT;
 
55
  YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(radius); 
 
56
  YGYOTO_WORKER_GETSET_OBJECT2(metric,Metric);
 
57
 
 
58
  /* INITCOORD */
 
59
  if ((iarg=kiargs[++k])>=0) { //initcoord
 
60
    iarg+=*rvset;
 
61
    if (yarg_nil(iarg)) {
 
62
      if ((*rvset)++) y_error(rmsg);
 
63
      long dims[]= {1, 8};
 
64
      double *coord = ypush_d(dims);
 
65
      (*ao)->getInitialCoord(coord);
 
66
    } else {
 
67
      long ptot=1, vtot=1;
 
68
      double *pos=ygeta_d(iarg, &ptot, 0);
 
69
      if (ptot<4) y_error("POS should have at least 4 elements");
 
70
      int n;
 
71
      double coord[8];
 
72
      for (n=0; n<4; ++n) coord[n]=pos[n];
 
73
 
 
74
      double *v=NULL, tdot0;
 
75
 
 
76
      if (yarg_number(piargs[0]+*rvset)) {
 
77
        if ((*paUsed)++) y_error(pmsg);
 
78
        v=ygeta_d(piargs[0]+*rvset, &vtot, 0);
 
79
        if (vtot!=3) y_error("V should have 3 elements");
 
80
        if (!(*ao)->metric())
 
81
          y_error("Please set metric before setting initial condition");
 
82
        tdot0=(*ao)->metric()->SysPrimeToTdot(pos, v);
 
83
        coord[4]=tdot0;
 
84
        for (n=0; n<3; ++n) coord[5+n]=v[n]*tdot0;
 
85
      } else if (ptot==8) {
 
86
        for (n=4; n<8; ++n) coord[n]=pos[n];
 
87
      } else if (ptot==7) {
 
88
        v=pos+4;
 
89
        tdot0=(*ao)->metric()->SysPrimeToTdot(pos, v);
 
90
        coord[4]=tdot0;
 
91
        for (n=0; n<3; ++n) coord[5+n]=v[n]*tdot0;
 
92
      } else y_error("Not enough information to set initial condition");
 
93
 
 
94
      (*ao)->setInitialCondition(coord);
 
95
    }
 
96
  }
 
97
 
 
98
  YGYOTO_WORKER_GETSET_OBJECT2(spectrum,Spectrum);
 
99
  YGYOTO_WORKER_GETSET_OBJECT2(opacity,Spectrum);
 
100
  YGYOTO_WORKER_GETSET_DOUBLE2_UNIT(delta);
 
101
  YGYOTO_WORKER_GETSET_LONG2(adaptive);
 
102
  YGYOTO_WORKER_GETSET_DOUBLE2(deltaMaxOverRadius);
 
103
  YGYOTO_WORKER_GETSET_DOUBLE2(deltaMaxOverDistance);
 
104
  YGYOTO_WORKER_GETSET_DOUBLE2(TMin);
 
105
  YGYOTO_WORKER_GETSET_DOUBLE2(TMax);
 
106
  YGYOTO_WORKER_GETSET_LONG2( maxiter );
 
107
  YGYOTO_WORKER_RUN( reset() );
 
108
  YGYOTO_WORKER_RUN( xFill(ygets_d(iarg+*rvset)) );
 
109
 
 
110
  YGYOTO_WORKER_CALL_GENERIC(Astrobj);
 
111
 
 
112
  k += YGYOTO_ASTROBJ_GENERIC_KW_N;
 
113
 
 
114
  // SPECIFIC GET KEYWORDS
 
115
  if ((iarg=kiargs[++k])>=0) { // get skypos
 
116
    if ((*rvset)++) y_error(rmsg);
 
117
    if ((*paUsed)++) y_error(pmsg);
 
118
 
 
119
    if (!yarg_Screen(iarg)) y_error("Expecting gyoto_Screen argument");
 
120
    SmartPointer<Screen> *screen = yget_Screen(iarg);
 
121
 
 
122
    int nel =(*ao)->get_nelements();
 
123
    long dims[] = {2, nel, 3};
 
124
    double * data=ypush_d(dims);
 
125
 
 
126
    (*ao)->getSkyPos(*screen, data, data+nel, data+2*nel);
 
127
  }
 
128
 
 
129
  if ((iarg=kiargs[++k])>=0) { // get_txyz
 
130
    if ((*rvset)++) y_error(rmsg);
 
131
    int nel =(*ao)->get_nelements();
 
132
      
 
133
    long dims[] = {2, nel, 4};
 
134
    double * data=ypush_d(dims);
 
135
 
 
136
    (*ao)->get_t(data);
 
137
    (*ao)->get_xyz(data+nel, data+2*nel, data+3*nel);
 
138
  }
 
139
 
 
140
  if ((iarg=kiargs[++k])>=0) { // get_prime
 
141
    if ((*rvset)++) y_error(rmsg);
 
142
    int nel =(*ao)->get_nelements();
 
143
      
 
144
    long dims[] = {2, nel, 3};
 
145
    double * data=ypush_d(dims);
 
146
 
 
147
    (*ao)->get_prime(data, data+nel, data+2*nel);
 
148
  }
 
149
 
 
150
  if ((iarg=kiargs[++k])>=0) { // get_coord
 
151
    if ((*rvset)++) y_error(rmsg);
 
152
    // Two cases : get_coord=1 or get_coord=dates. 
 
153
    // We will recognize the first case if the parameter is the
 
154
    // integer scalar 1, any other type or value is a date
 
155
    // specification.
 
156
    int argt = yarg_number(iarg);
 
157
    long ntot[1] = {1};
 
158
    long dims[Y_DIMSIZE];
 
159
    double* dates = NULL;
 
160
    double dummy[]= {1.};
 
161
    if (yarg_nil(iarg)) {
 
162
      // void is same as 1 in this context
 
163
      argt=1;
 
164
      *ntot=1;
 
165
      dims[0]=0;
 
166
      dates=dummy;
 
167
    } else dates = ygeta_d(iarg, ntot, dims);
 
168
 
 
169
    if (debug())
 
170
      cerr << "DEBUG: gyoto_StarTrace(get_coord=array("
 
171
           << (argt==1?"long":"double")
 
172
           << ", " << *ntot 
 
173
           << ") (rank=" << dims[0]
 
174
           << ", dates[0]=" << dates[0] << ")" << endl;
 
175
 
 
176
    if (dims[0] == 0 && argt == 1 && dates[0] == 1) {
 
177
      if(debug())
 
178
        cerr << "DEBUG: retrieving all already computed coordinates" << endl;
 
179
      int nel =(*ao)->get_nelements();
 
180
      long ddims[] = {2, nel, 8};
 
181
      double * data = ypush_d(ddims);
 
182
 
 
183
      (*ao)->getCoord(data, data+nel, data+2*nel, data+3*nel);
 
184
      (*ao)->get_dot(data+4*nel, data+5*nel, data+6*nel, data+7*nel);
 
185
    } else {
 
186
      if(debug())
 
187
        cerr << "DEBUG: retrieving coordinates for specified date(s)"<< endl;
 
188
      if (dims[0] > Y_DIMSIZE-2)
 
189
        y_errorn("gyoto_StarTrace(get_coord=dates): DATES must have at most %ld "
 
190
                 "dimensions", Y_DIMSIZE-2);
 
191
      dims[0] += 1;
 
192
      dims[dims[0]] = 7;
 
193
      double * data = ypush_d(dims);
 
194
      size_t nel = *ntot;
 
195
      (*ao) -> getCoord(dates, nel,
 
196
                        data, data+ nel,data+2* nel,data+3* nel,
 
197
                        data+4* nel, data+5* nel, data+6* nel);
 
198
    }
 
199
  }
 
200
  if ((iarg=kiargs[++k])>=0) { // get_cartesian
 
201
    if ((*rvset)++) y_error(rmsg);
 
202
    long ntot[1] = {1};
 
203
    long dims[Y_DIMSIZE];
 
204
    double* dates = ygeta_d(iarg, ntot, dims);
 
205
    if(debug())
 
206
      cerr << "DEBUG: gyoto_StarTrace(get_cartesian=dates)"<< endl;
 
207
    if (dims[0] > Y_DIMSIZE-2)
 
208
      y_errorn("gyoto_StarTrace(get_cartesian=dates): DATES must have at most %ld "
 
209
               "dimensions", Y_DIMSIZE-2);
 
210
    dims[0] += 1;
 
211
    dims[dims[0]] = 6;
 
212
    double * data = ypush_d(dims);
 
213
    size_t nel = *ntot;
 
214
    (*ao) -> getCartesian(dates, nel,
 
215
                          data,         data +   nel, data + 2*nel,
 
216
                          data + 3*nel, data + 4*nel, data + 5*nel);
 
217
    
 
218
  }
 
219
 
 
220
  if ((iarg=kiargs[++k])>=0) { // star
 
221
    if ((*rvset)++) y_error(rmsg);
 
222
    *ypush_Astrobj() = new Star (**ao);
 
223
  }
 
224
  
 
225
}
 
226
 
 
227
 
 
228
 
 
229
extern "C" {
 
230
  void Y__gyoto_StarTrace_register_as_Astrobj(){
 
231
    ygyoto_Astrobj_register("StarTrace",&ygyoto_StarTrace_eval);
 
232
  }
 
233
 
 
234
  // Generic constructor/accessor
 
235
  void
 
236
  Y_gyoto_StarTrace(int argc)
 
237
  {
 
238
    GYOTO_DEBUG << endl;
 
239
    YGYOTO_CONSTRUCTOR_INIT2(Astrobj, Astrobj::Generic, StarTrace, astrobj);
 
240
    if ((*ao)->kind().compare("StarTrace"))
 
241
      y_error("Expecting Astrobj of kind StarTrace");
 
242
    ygyoto_StarTrace_eval(ao, argc);
 
243
  }
 
244
 
 
245
}