source: CIVL/examples/omp/shtns/shtns_numpy.i

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 24.2 KB
RevLine 
[2131108]1/*
2 * Copyright (c) 2010-2013 Centre National de la Recherche Scientifique.
3 * written by Nathanael Schaeffer (CNRS, ISTerre, Grenoble, France).
4 *
5 * nathanael.schaeffer@ujf-grenoble.fr
6 *
7 * This software is governed by the CeCILL license under French law and
8 * abiding by the rules of distribution of free software. You can use,
9 * modify and/or redistribute the software under the terms of the CeCILL
10 * license as circulated by CEA, CNRS and INRIA at the following URL
11 * "http://www.cecill.info".
12 *
13 * The fact that you are presently reading this means that you have had
14 * knowledge of the CeCILL license and that you accept its terms.
15 *
16 */
17
18/* shtns_numpy.i : SWIG interface to Python using NumPy */
19
20/* TODO and known problems :
21 * - alignement on 16 bytes of NumPy arrays is not guaranteed. It should however work on 64bit systems or on modern 32bit systems.
22 * - you may have to adjust the path below to include the header file "arrayobject.h" from the NumPy package.
23 */
24
25%module (docstring="Python/NumPy interface to the SHTns spherical harmonic transform library") shtns
26
27%init{
28 import_array(); // required by NumPy
29}
30
31%pythoncode{
32 import numpy as np
33}
34
35%{
36
37#include <numpy/arrayobject.h>
38#include "sht_private.h"
39
40// variables used for exception handling.
41static int shtns_error = 0;
42static char* shtns_err_msg;
43static char msg_buffer[128];
44static char msg_grid_err[] = "Grid not set. Call .set_grid() mehtod.";
45static char msg_numpy_arr[] = "Numpy array expected.";
46static char msg_rot_err[] = "truncation must be triangular (lmax=mmax, mres=1)";
47
48static void throw_exception(int error, int iarg, char* msg)
49{
50 shtns_error = error;
51 shtns_err_msg = msg;
52 if (iarg > 0) {
53 sprintf(msg_buffer, "arg #%d : %.100s", iarg, msg);
54 shtns_err_msg = msg_buffer;
55 }
56}
57
58static int check_spatial(int i, PyObject *a, int size) {
59 if (size == 0) {
60 throw_exception(SWIG_RuntimeError,0,msg_grid_err); return 0;
61 }
62 if (!PyArray_Check(a)) {
63 throw_exception(SWIG_TypeError,i,msg_numpy_arr); return 0;
64 }
65 if (PyArray_TYPE(a) != PyArray_DOUBLE) {
66 throw_exception(SWIG_TypeError,i,"spatial array must consist of float."); return 0;
67 }
68 if (!PyArray_ISCONTIGUOUS(a)) {
69 throw_exception(SWIG_RuntimeError,i,"spatial array not contiguous. Use 'b=a.copy()' to copy a to a contiguous array b."); return 0;
70 }
71 if (PyArray_SIZE(a) != size) {
72 throw_exception(SWIG_RuntimeError,i,"spatial array has wrong size"); return 0;
73 }
74 return 1;
75}
76
77static int check_spectral(int i, PyObject *a, int size) {
78 if (!PyArray_Check(a)) {
79 throw_exception(SWIG_RuntimeError,i,msg_numpy_arr); return 0;
80 }
81 if (PyArray_TYPE(a) != PyArray_CDOUBLE) {
82 throw_exception(SWIG_RuntimeError,i,"spectral array must consist of complex float. Create with: 'sh.spec_array()'"); return 0;
83 }
84 if (!PyArray_ISCONTIGUOUS(a)) {
85 throw_exception(SWIG_RuntimeError,i,"spactral array not contiguous. Use .copy() to copy to a contiguous array."); return 0;
86 }
87 if (PyArray_SIZE(a) != size) {
88 throw_exception(SWIG_RuntimeError,i,"spectral array has wrong size"); return 0;
89 }
90 return 1;
91}
92
93inline PyObject* SpecArray_New(int size) {
94 npy_intp dims = size;
95 npy_intp strides = sizeof(cplx);
96 return PyArray_New(&PyArray_Type, 1, &dims, PyArray_CDOUBLE, &strides, NULL, strides, 0, NULL);
97}
98
99inline PyObject* SpatArray_New(int size) {
100 npy_intp dims = size;
101 npy_intp strides = sizeof(double);
102 return PyArray_New(&PyArray_Type, 1, &dims, PyArray_DOUBLE, &strides, NULL, strides, 0, NULL);
103}
104
105%}
106
107// main object is renamed to sht.
108%rename("sht") shtns_info;
109%ignore SHT_NATIVE_LAYOUT;
110%ignore nlat_2;
111%ignore lmidx;
112%ignore li;
113%ignore mi;
114%ignore ct;
115%ignore st;
116
117%rename(print_version) shtns_print_version;
118%rename(set_verbosity) shtns_verbose;
119
120%feature("autodoc");
121%include "shtns.h"
122%include "exception.i"
123
124
125%extend shtns_info {
126 %exception {
127 shtns_error = 0; // clear exception
128 $function
129 if (shtns_error) { // test for exception
130 SWIG_exception(shtns_error, shtns_err_msg); return NULL;
131 }
132 }
133
134 %pythonappend shtns_info %{
135 ## array giving the degree of spherical harmonic coefficients.
136 self.l = np.zeros(self.nlm, dtype=np.int32)
137 ## array giving the order of spherical harmonic coefficients.
138 self.m = np.zeros(self.nlm, dtype=np.int32)
139 for mloop in range(0, self.mmax*self.mres+1, self.mres):
140 for lloop in range(mloop, self.lmax+1):
141 ii = self.idx(lloop,mloop)
142 self.m[ii] = mloop
143 self.l[ii] = lloop
144 self.m.flags.writeable = False # prevent writing in m and l arrays
145 self.l.flags.writeable = False
146 %}
147 %feature("kwargs") shtns_info;
148 shtns_info(int lmax, int mmax=-1, int mres=1, int norm=sht_orthonormal, int nthreads=0) { // default arguments : mmax, mres and norm
149 if (lmax < 2) {
150 throw_exception(SWIG_ValueError,1,"lmax < 2 not allowed"); return NULL;
151 }
152 if (mres <= 0) {
153 throw_exception(SWIG_ValueError,3,"mres <= 0 invalid"); return NULL;
154 }
155 if (mmax < 0) mmax = lmax/mres; // default mmax
156 if (mmax*mres > lmax) {
157 throw_exception(SWIG_ValueError,1,"lmax < mmax*mres invalid"); return NULL;
158 }
159 shtns_use_threads(nthreads); // use nthreads openmp threads if available (0 means auto)
160 return shtns_create(lmax, mmax, mres, norm);
161 }
162
163 ~shtns_info() {
164 shtns_destroy($self); // free memory.
165 }
166
167 %pythonappend set_grid %{
168 ## array giving the cosine of the colatitude for the grid.
169 self.cos_theta = self.__ct()
170 self.cos_theta.flags.writeable = False
171 ## shape of a spatial array for the grid (tuple of 2 values).
172 self.spat_shape = tuple(self.__spat_shape())
173 %}
174 %apply int *OUTPUT { int *nlat_out };
175 %apply int *OUTPUT { int *nphi_out };
176 %feature("kwargs") set_grid;
177 void set_grid(int nlat=0, int nphi=0, int flags=sht_quick_init, double polar_opt=1.0e-8, int nl_order=1, int *nlat_out, int *nphi_out) { // default arguments
178 if (nlat != 0) {
179 if (nlat <= $self->lmax) { // nlat too small
180 throw_exception(SWIG_ValueError,1,"nlat <= lmax"); return;
181 }
182 if (nlat & 1) { // nlat must be even
183 throw_exception(SWIG_ValueError,1,"nlat must be even"); return;
184 }
185 }
186 if ((nphi != 0) && (nphi <= $self->mmax *2)) { // nphi too small
187 throw_exception(SWIG_ValueError,2,"nphi <= 2*mmax"); return;
188 }
189 if (!(flags & SHT_THETA_CONTIGUOUS)) flags |= SHT_PHI_CONTIGUOUS; // default to SHT_PHI_CONTIGUOUS.
190 *nlat_out = nlat; *nphi_out = nphi;
191 shtns_set_grid_auto($self, flags, polar_opt, nl_order, nlat_out, nphi_out);
192 }
193
194 void print_info() {
195 shtns_print_cfg($self);
196 }
197 double sh00_1() {
198 return sh00_1($self);
199 }
200 double sh10_ct() {
201 return sh10_ct($self);
202 }
203 double sh11_st() {
204 return sh11_st($self);
205 }
206 double shlm_e1(unsigned l, unsigned m) {
207 return shlm_e1($self, l, m);
208 }
209
210 /* returns useful data */
211 PyObject* __ct() { // grid must have been initialized.
212 PyObject *obj;
213 double *ct;
214 int i;
215 if ($self->nlat == 0) { // no grid
216 throw_exception(SWIG_RuntimeError,0,msg_grid_err);
217 return NULL;
218 }
219 obj = SpatArray_New($self->nlat);
220 ct = (double*) PyArray_DATA(obj);
221 for (i=0; i<$self->nlat; i++) ct[i] = $self->ct[i]; // copy
222 return obj;
223 }
224 PyObject* gauss_wts() { // gauss grid must have been initialized.
225 PyObject *obj;
226 if ($self->nlat == 0) { // no grid
227 throw_exception(SWIG_RuntimeError,0,msg_grid_err);
228 return NULL;
229 }
230 if ($self->wg == NULL) {
231 throw_exception(SWIG_RuntimeError,0,"not a gauss grid");
232 return NULL;
233 }
234 obj = SpatArray_New($self->nlat_2);
235 shtns_gauss_wts($self, PyArray_DATA(obj));
236 return obj;
237 }
238 PyObject* mul_ct_matrix() {
239 PyObject *mx = SpatArray_New(2*$self->nlm);
240 mul_ct_matrix($self, PyArray_DATA(mx));
241 return mx;
242 }
243 PyObject* st_dt_matrix() {
244 PyObject *mx = SpatArray_New(2*$self->nlm);
245 st_dt_matrix($self, PyArray_DATA(mx));
246 return mx;
247 }
248
249 %apply int *OUTPUT { int *dim0 };
250 %apply int *OUTPUT { int *dim1 };
251 void __spat_shape(int *dim0, int *dim1) {
252 *dim0 = $self->nphi; *dim1 = $self->nlat;
253 if ($self->fftc_mode == 1) { // phi-contiguous
254 *dim0 = $self->nlat; *dim1 = $self->nphi;
255 }
256 }
257
258 %pythoncode %{
259 def spec_array(self, im=-1):
260 """return a numpy array of spherical harmonic coefficients (complex). Adress coefficients with index sh.idx(l,m)
261 if optional argument im is given, the spectral array is restricted to order im*mres."""
262 if im<0:
263 return np.zeros(self.nlm, dtype=complex)
264 else:
265 return np.zeros(self.lmax + 1 - im*self.mres, dtype=complex)
266
267 def spat_array(self):
268 """return a numpy array of 2D spatial field."""
269 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
270 return np.zeros(self.spat_shape)
271 %}
272
273 // returns the index in a spectral array of (l,m) coefficient.
274 int idx(unsigned l, unsigned m) {
275 if (l > $self->lmax) {
276 throw_exception(SWIG_ValueError,1,"l invalid"); return 0;
277 }
278 if ( (m > l) || (m > $self->mmax * $self->mres) || (m % $self->mres != 0) ) {
279 throw_exception(SWIG_ValueError,2,"m invalid"); return 0;
280 }
281 return LM($self, l, m);
282 }
283
284 /* scalar transforms */
285 void spat_to_SH(PyObject *Vr, PyObject *Qlm) {
286 if (check_spatial(1,Vr, $self->nspat) && check_spectral(2,Qlm, $self->nlm))
287 spat_to_SH($self, PyArray_DATA(Vr), PyArray_DATA(Qlm));
288 }
289 void SH_to_spat(PyObject *Qlm, PyObject *Vr) {
290 if (check_spatial(2,Vr, $self->nspat) && check_spectral(1,Qlm, $self->nlm))
291 SH_to_spat($self, PyArray_DATA(Qlm), PyArray_DATA(Vr));
292 }
293 /* complex transforms */
294 void spat_cplx_to_SH(PyObject *z, PyObject *alm) {
295 int n = $self->lmax + 1;
296 if (check_spectral(1,z, $self->nspat) && check_spectral(2,alm, n*n))
297 spat_cplx_to_SH($self, PyArray_DATA(z), PyArray_DATA(alm));
298 }
299 void SH_to_spat_cplx(PyObject *alm, PyObject *z) {
300 int n = $self->lmax + 1;
301 if (check_spectral(2,z, $self->nspat) && check_spectral(1,alm, n*n))
302 SH_to_spat_cplx($self, PyArray_DATA(alm), PyArray_DATA(z));
303 }
304 /*void SH_to_spat_grad(PyObject *alm, PyObject *gt, PyObject *gp) {
305 if (check_spatial(3,gp, $self->nspat) && check_spatial(2,gt, $self->nspat) && check_spectral(1,alm, $self->nlm))
306 SH_to_spat_grad($self, PyArray_DATA(alm), PyArray_DATA(gt), PyArray_DATA(gp));
307 }*/
308
309 /* 2D vectors */
310 void spat_to_SHsphtor(PyObject *Vt, PyObject *Vp, PyObject *Slm, PyObject *Tlm) {
311 if (check_spatial(1,Vt, $self->nspat) && check_spatial(2,Vp, $self->nspat) && check_spectral(3,Slm, $self->nlm) && check_spectral(4,Tlm, $self->nlm))
312 spat_to_SHsphtor($self, PyArray_DATA(Vt), PyArray_DATA(Vp), PyArray_DATA(Slm), PyArray_DATA(Tlm));
313 }
314 void SHsphtor_to_spat(PyObject *Slm, PyObject *Tlm, PyObject *Vt, PyObject *Vp) {
315 if (check_spatial(3,Vt, $self->nspat) && check_spatial(4,Vp, $self->nspat) && check_spectral(1,Slm, $self->nlm) && check_spectral(2,Tlm, $self->nlm))
316 SHsphtor_to_spat($self, PyArray_DATA(Slm), PyArray_DATA(Tlm), PyArray_DATA(Vt), PyArray_DATA(Vp));
317 }
318 void SHsph_to_spat(PyObject *Slm, PyObject *Vt, PyObject *Vp) {
319 if (check_spatial(2,Vt, $self->nspat) && check_spatial(3,Vp, $self->nspat) && check_spectral(1,Slm, $self->nlm))
320 SHsph_to_spat($self, PyArray_DATA(Slm), PyArray_DATA(Vt), PyArray_DATA(Vp));
321 }
322 void SHtor_to_spat(PyObject *Tlm, PyObject *Vt, PyObject *Vp) {
323 if (check_spatial(2,Vt, $self->nspat) && check_spatial(3,Vp, $self->nspat) && check_spectral(1,Tlm, $self->nlm))
324 SHtor_to_spat($self, PyArray_DATA(Tlm), PyArray_DATA(Vt), PyArray_DATA(Vp));
325 }
326
327 /* 3D vectors */
328 void spat_to_SHqst(PyObject *Vr, PyObject *Vt, PyObject *Vp, PyObject *Qlm, PyObject *Slm, PyObject *Tlm) {
329 if (check_spatial(1,Vr, $self->nspat) && check_spatial(2,Vt, $self->nspat) && check_spatial(3,Vp, $self->nspat)
330 && check_spectral(4,Qlm, $self->nlm) && check_spectral(5,Slm, $self->nlm) && check_spectral(6,Tlm, $self->nlm))
331 spat_to_SHqst($self, PyArray_DATA(Vr), PyArray_DATA(Vt), PyArray_DATA(Vp), PyArray_DATA(Qlm), PyArray_DATA(Slm), PyArray_DATA(Tlm));
332 }
333 void SHqst_to_spat(PyObject *Qlm, PyObject *Slm, PyObject *Tlm, PyObject *Vr, PyObject *Vt, PyObject *Vp) {
334 if (check_spatial(4,Vr, $self->nspat) && check_spatial(5,Vt, $self->nspat) && check_spatial(6,Vp, $self->nspat)
335 && check_spectral(1,Qlm, $self->nlm) && check_spectral(2,Slm, $self->nlm) && check_spectral(3,Tlm, $self->nlm))
336 SHqst_to_spat($self, PyArray_DATA(Qlm), PyArray_DATA(Slm), PyArray_DATA(Tlm), PyArray_DATA(Vr), PyArray_DATA(Vt), PyArray_DATA(Vp));
337 }
338
339 %pythoncode %{
340 def synth(self,*arg):
341 """
342 spectral to spatial transform, for scalar or vector data.
343 v = synth(qlm) : compute the spatial representation of the scalar qlm
344 vtheta,vphi = synth(slm,tlm) : compute the 2D spatial vector from its spectral spheroidal/toroidal scalars (slm,tlm)
345 vr,vtheta,vphi = synth(qlm,slm,tlm) : compute the 3D spatial vector from its spectral radial/spheroidal/toroidal scalars (qlm,slm,tlm)
346 """
347 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
348 n = len(arg)
349 if (n>3) or (n<1): raise RuntimeError("1,2 or 3 arguments required.")
350 q = list(arg)
351 for i in range(0,n):
352 if q[i].size != self.nlm: raise RuntimeError("spectral array has wrong size.")
353 if q[i].dtype.num != np.dtype('complex128').num: raise RuntimeError("spectral array should be dtype=complex.")
354 if q[i].flags.contiguous == False: q[i] = q[i].copy() # contiguous array required.
355 if n==1: #scalar transform
356 vr = np.empty(self.spat_shape)
357 self.SH_to_spat(q[0],vr)
358 return vr
359 elif n==2: # 2D vector transform
360 vt = np.empty(self.spat_shape) # v_theta
361 vp = np.empty(self.spat_shape) # v_phi
362 self.SHsphtor_to_spat(q[0],q[1],vt,vp)
363 return vt,vp
364 else: # 3D vector transform
365 vr = np.empty(self.spat_shape) # v_r
366 vt = np.empty(self.spat_shape) # v_theta
367 vp = np.empty(self.spat_shape) # v_phi
368 self.SHqst_to_spat(q[0],q[1],q[2],vr,vt,vp)
369 return vr,vt,vp
370
371 def analys(self,*arg):
372 """
373 spatial to spectral transform, for scalar or vector data.
374 qlm = analys(q) : compute the spherical harmonic representation of the scalar q
375 slm,tlm = analys(vtheta,vphi) : compute the spectral spheroidal/toroidal scalars (slm,tlm) from 2D vector components (vtheta, vphi)
376 qlm,slm,tlm = synth(vr,vtheta,vphi) : compute the spectral radial/spheroidal/toroidal scalars (qlm,slm,tlm) from 3D vector components (vr,vtheta,vphi)
377 """
378 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
379 if abs(self.cos_theta[0]) == 1: raise RuntimeError("Analysis not allowed with sht_reg_poles grid.")
380 n = len(arg)
381 if (n>3) or (n<1): raise RuntimeError("1,2 or 3 arguments required.")
382 v = list(arg)
383 for i in range(0,n):
384 if v[i].shape != self.spat_shape: raise RuntimeError("spatial array has wrong shape.")
385 if v[i].dtype.num != np.dtype('float64').num: raise RuntimeError("spatial array should be dtype=float64.")
386 if v[i].flags.contiguous == False: v[i] = v[i].copy() # contiguous array required.
387 if n==1:
388 q = np.empty(self.nlm, dtype=complex)
389 self.spat_to_SH(v[0],q)
390 return q
391 elif n==2:
392 s = np.empty(self.nlm, dtype=complex)
393 t = np.empty(self.nlm, dtype=complex)
394 self.spat_to_SHsphtor(v[0],v[1],s,t)
395 return s,t
396 else:
397 q = np.empty(self.nlm, dtype=complex)
398 s = np.empty(self.nlm, dtype=complex)
399 t = np.empty(self.nlm, dtype=complex)
400 self.spat_to_SHqst(v[0],v[1],v[2],q,s,t)
401 return q,s,t
402
403 def synth_grad(self,slm):
404 """(vtheta,vphi) = synth_grad(sht self, slm) : compute the spatial representation of the gradient of slm"""
405 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
406 if slm.size != self.nlm: raise RuntimeError("spectral array has wrong size.")
407 if slm.dtype.num != np.dtype('complex128').num: raise RuntimeError("spectral array should be dtype=complex.")
408 if slm.flags.contiguous == False: slm = slm.copy() # contiguous array required.
409 vt = np.empty(self.spat_shape)
410 vp = np.empty(self.spat_shape)
411 self.SHsph_to_spat(slm,vt,vp)
412 return vt,vp
413
414 def synth_cplx(self,alm):
415 """
416 spectral to spatial transform, for complex valued scalar data.
417 z = synth(alm) : compute the spatial representation of the scalar alm
418 """
419 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
420 if self.lmax != self.mmax: raise RuntimeError("complex SH requires lmax=mmax and mres=1.")
421 if alm.size != (self.lmax+1)**2: raise RuntimeError("spectral array has wrong size.")
422 if alm.dtype.num != np.dtype('complex128').num: raise RuntimeError("spectral array should be dtype=complex.")
423 if alm.flags.contiguous == False: alm = alm.copy() # contiguous array required.
424 z = np.empty(self.spat_shape, dtype=complex)
425 self.SH_to_spat_cplx(alm,z)
426 return z
427
428 def analys_cplx(self,z):
429 """
430 spatial to spectral transform, for complex valued scalar data.
431 alm = analys(z) : compute the spherical harmonic representation of the complex scalar z
432 """
433 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
434 if self.lmax != self.mmax: raise RuntimeError("complex SH requires lmax=mmax and mres=1.")
435 if z.shape != self.spat_shape: raise RuntimeError("spatial array has wrong shape.")
436 if z.dtype.num != np.dtype('complex128').num: raise RuntimeError("spatial array should be dtype=complex128.")
437 if z.flags.contiguous == False: z = z.copy() # contiguous array required.
438 alm = np.empty((self.lmax+1)**2, dtype=complex)
439 self.spat_cplx_to_SH(z,alm)
440 return alm
441
442 def zidx(self, l,m):
443 """
444 zidx(sht self, int l, int m) -> int : compute the index l*(l+1)+m in a complex spherical harmonic expansion
445 """
446 l = np.asarray(l)
447 m = np.asarray(m)
448 if (l>self.lmax).any() or (abs(m)>l).any() : raise RuntimeError("invalid range for l,m")
449 return l*(l+1)+m
450
451 def zlm(self, idx):
452 """
453 zlm(sht self, int idx) -> (int,int) : returns the l and m corresponding to the given index in complex spherical harmonic expansion
454 """
455 idx = np.asarray(idx)
456 if (idx >= (self.lmax+1)**2).any() or (idx < 0).any() : raise RuntimeError("invalid range for l,m")
457 l = np.sqrt(idx).astype(int)
458 m = idx - l*(l+1)
459 return l,m
460
461 def spec_array_cplx(self):
462 """return a numpy array that can hold the spectral representation of a complex scalar spatial field."""
463 return np.zeros((self.lmax+1)**2, dtype=complex)
464
465 def spat_array_cplx(self):
466 """return a numpy array of 2D complex spatial field."""
467 if self.nlat == 0: raise RuntimeError("Grid not set. Call .set_grid() mehtod.")
468 return np.zeros(self.spat_shape, dtype=complex128)
469 %}
470
471 /* local evaluations */
472 double SH_to_point(PyObject *Qlm, double cost, double phi) {
473 if (check_spectral(1,Qlm, $self->nlm)) return SH_to_point($self, PyArray_DATA(Qlm), cost, phi);
474 return 0.0;
475 }
476 %apply double *OUTPUT { double *vr };
477 %apply double *OUTPUT { double *vt };
478 %apply double *OUTPUT { double *vp };
479 void SH_to_grad_point(PyObject *DrSlm, PyObject *Slm, double cost, double phi, double *vr, double *vt, double *vp) {
480 if (check_spectral(1,DrSlm, $self->nlm) && check_spectral(2,Slm, $self->nlm))
481 SH_to_grad_point($self, PyArray_DATA(DrSlm), PyArray_DATA(Slm), cost, phi, vr, vt, vp);
482 }
483 void SHqst_to_point(PyObject *Qlm, PyObject *Slm, PyObject *Tlm,
484 double cost, double phi, double *vr, double *vt, double *vp) {
485 if (check_spectral(1,Qlm, $self->nlm) && check_spectral(2,Slm, $self->nlm) && check_spectral(3,Tlm, $self->nlm))
486 SHqst_to_point($self, PyArray_DATA(Qlm), PyArray_DATA(Slm), PyArray_DATA(Tlm), cost, phi, vr, vt, vp);
487 }
488 %clear double *vr;
489 %clear double *vt;
490 %clear double *vp;
491
492 /* rotation of SH representations (experimental) */
493 PyObject* Zrotate(PyObject *Qlm, double alpha) {
494 if (check_spectral(1,Qlm, $self->nlm)) {
495 PyObject *Rlm = SpecArray_New($self->nlm);
496 SH_Zrotate($self, PyArray_DATA(Qlm), alpha, PyArray_DATA(Rlm));
497 return Rlm;
498 }
499 return NULL;
500 }
501 PyObject* Yrotate(PyObject *Qlm, double alpha) {
502 if (($self->mres != 1)||($self->mmax != $self->lmax)) {
503 throw_exception(SWIG_RuntimeError,0,msg_rot_err); return NULL;
504 }
505 if (check_spectral(1,Qlm, $self->nlm)) {
506 PyObject *Rlm = SpecArray_New($self->nlm);
507 SH_Yrotate($self, PyArray_DATA(Qlm), alpha, PyArray_DATA(Rlm));
508 return Rlm;
509 }
510 return NULL;
511 }
512 PyObject* Yrotate90(PyObject *Qlm) {
513 if (($self->mres != 1)||($self->mmax != $self->lmax)) {
514 throw_exception(SWIG_RuntimeError,0,msg_rot_err); return NULL;
515 }
516 if (check_spectral(1,Qlm, $self->nlm)) {
517 PyObject *Rlm = SpecArray_New($self->nlm);
518 SH_Yrotate90($self, PyArray_DATA(Qlm), PyArray_DATA(Rlm));
519 return Rlm;
520 }
521 return NULL;
522 }
523 PyObject* Xrotate90(PyObject *Qlm) {
524 if (($self->mres != 1)||($self->mmax != $self->lmax)) {
525 throw_exception(SWIG_RuntimeError,0,msg_rot_err); return NULL;
526 }
527 if (check_spectral(1,Qlm, $self->nlm)) {
528 PyObject *Rlm = SpecArray_New($self->nlm);
529 SH_Xrotate90($self, PyArray_DATA(Qlm), PyArray_DATA(Rlm));
530 return Rlm;
531 }
532 return NULL;
533 }
534
535 /* multiplication by l+1 l-1 matrix (mul_ct_matrix or st_dt_matrix) */
536 PyObject* SH_mul_mx(PyObject *mx, PyObject *Qlm) {
537 if (check_spectral(2,Qlm, $self->nlm) && check_spatial(1, mx, 2* $self->nlm)) {
538 PyObject *Rlm = SpecArray_New($self->nlm);
539 SH_mul_mx($self, PyArray_DATA(mx), PyArray_DATA(Qlm), PyArray_DATA(Rlm));
540 return Rlm;
541 }
542 return NULL;
543 }
544
545 /* Legendre transforms (no fft) at given order m */
546 void spat_to_SH_m(PyObject *Vr, PyObject *Qlm, PyObject *im) {
547 int im_ = PyLong_AsLong(im); int ltr = $self->lmax;
548 if ((im_ >= 0) && check_spectral(1,Vr, $self->nlat) && check_spectral(2,Qlm, ltr+1 - im_*$self->mres))
549 spat_to_SH_ml($self, im_, PyArray_DATA(Vr), PyArray_DATA(Qlm), ltr);
550 }
551 void SH_to_spat_m(PyObject *Qlm, PyObject *Vr, PyObject *im) {
552 int im_ = PyLong_AsLong(im); int ltr = $self->lmax;
553 if ((im_ >= 0) && check_spectral(2,Vr, $self->nlat) && check_spectral(1,Qlm, ltr+1 - im_*$self->mres))
554 SH_to_spat_ml($self, im_, PyArray_DATA(Qlm), PyArray_DATA(Vr), ltr);
555 }
556 void spat_to_SHsphtor_m(PyObject *Vt, PyObject *Vp, PyObject *Slm, PyObject *Tlm, PyObject *im) {
557 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
558 if ((im_ >= 0) && check_spectral(1,Vt, $self->nlat) && check_spectral(2,Vp, $self->nlat) && check_spectral(3,Slm, nelem) && check_spectral(4,Tlm, nelem))
559 spat_to_SHsphtor_ml($self, im_, PyArray_DATA(Vt), PyArray_DATA(Vp), PyArray_DATA(Slm), PyArray_DATA(Tlm), ltr);
560 }
561 void SHsphtor_to_spat_m(PyObject *Slm, PyObject *Tlm, PyObject *Vt, PyObject *Vp, PyObject *im) {
562 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
563 if ((im_ >= 0) && check_spectral(3,Vt, $self->nlat) && check_spectral(4,Vp, $self->nlat) && check_spectral(1,Slm, nelem) && check_spectral(2,Tlm, nelem))
564 SHsphtor_to_spat_ml($self, im_, PyArray_DATA(Slm), PyArray_DATA(Tlm), PyArray_DATA(Vt), PyArray_DATA(Vp), ltr);
565 }
566 void SHsph_to_spat_m(PyObject *Slm, PyObject *Vt, PyObject *Vp, PyObject *im) {
567 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
568 if ((im_ >= 0) && check_spectral(2,Vt, $self->nlat) && check_spectral(3,Vp, $self->nlat) && check_spectral(1,Slm, nelem))
569 SHsph_to_spat_ml($self, im_, PyArray_DATA(Slm), PyArray_DATA(Vt), PyArray_DATA(Vp), ltr);
570 }
571 void SHtor_to_spat_m(PyObject *Tlm, PyObject *Vt, PyObject *Vp, PyObject *im) {
572 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
573 if ((im_ >= 0) && check_spectral(2,Vt, $self->nlat) && check_spectral(3,Vp, $self->nlat) && check_spectral(1,Tlm, nelem))
574 SHtor_to_spat_ml($self, im_, PyArray_DATA(Tlm), PyArray_DATA(Vt), PyArray_DATA(Vp), ltr);
575 }
576 void spat_to_SHqst_m(PyObject *Vr, PyObject *Vt, PyObject *Vp, PyObject *Qlm, PyObject *Slm, PyObject *Tlm, PyObject *im) {
577 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
578 if ((im_ >= 0) && check_spectral(1,Vr, $self->nlat) && check_spectral(2,Vt, $self->nlat) && check_spectral(3,Vp, $self->nlat)
579 && check_spectral(4,Qlm, nelem) && check_spectral(5,Slm, nelem) && check_spectral(6,Tlm, nelem))
580 spat_to_SHqst_ml($self, im_, PyArray_DATA(Vr), PyArray_DATA(Vt), PyArray_DATA(Vp), PyArray_DATA(Qlm), PyArray_DATA(Slm), PyArray_DATA(Tlm), ltr);
581 }
582 void SHqst_to_spat_m(PyObject *Qlm, PyObject *Slm, PyObject *Tlm, PyObject *Vr, PyObject *Vt, PyObject *Vp, PyObject *im) {
583 int im_ = PyLong_AsLong(im); int ltr = $self->lmax; int nelem = ltr+1 - im_*$self->mres;
584 if ((im_ >= 0) && check_spectral(4,Vr, $self->nlat) && check_spectral(5,Vt, $self->nlat) && check_spectral(6,Vp, $self->nlat)
585 && check_spectral(1,Qlm, nelem) && check_spectral(2,Slm, nelem) && check_spectral(3,Tlm, nelem))
586 SHqst_to_spat_ml($self, im_, PyArray_DATA(Qlm), PyArray_DATA(Slm), PyArray_DATA(Tlm), PyArray_DATA(Vr), PyArray_DATA(Vt), PyArray_DATA(Vp), ltr);
587 }
588
589};
Note: See TracBrowser for help on using the repository browser.