← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/specfun_wrappers.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/specfun_wrappers.c
 
1
/* This file is a collection of wrappers around the
 
2
 *  Specfun Fortran library of functions 
 
3
 */
 
4
 
 
5
#include "specfun_wrappers.h"
 
6
#include <stdio.h>
 
7
 
 
8
#define CADDR(z) (double *)(&((z).real)), (double*)(&((z).imag))
 
9
#define F2C_CST(z) (double *)&((z)->real), (double *)&((z)->imag)
 
10
 
 
11
#if defined(NO_APPEND_FORTRAN)
 
12
#if defined(UPPERCASE_FORTRAN)
 
13
#define F_FUNC(f,F) F
 
14
#else
 
15
#define F_FUNC(f,F) f
 
16
#endif
 
17
#else
 
18
#if defined(UPPERCASE_FORTRAN)
 
19
#define F_FUNC(f,F) F##_
 
20
#else
 
21
#define F_FUNC(f,F) f##_
 
22
#endif
 
23
#endif
 
24
 
 
25
extern double psi(double);
 
26
extern double struve(double, double);
 
27
 
 
28
extern void F_FUNC(cgama,CGAMA)(double*,double*,int*,double*,double*);
 
29
extern void F_FUNC(cpsi,CPSI)(double*,double*,double*,double*);
 
30
extern void F_FUNC(hygfz,HYGFZ)(double*,double*,double*,Py_complex*,Py_complex*);
 
31
extern void F_FUNC(cchg,CCHG)(double*,double*,Py_complex*,Py_complex*);
 
32
extern void F_FUNC(chgu,CHGU)(double*,double*,double*,double*,int*);
 
33
extern void F_FUNC(itairy,ITAIRY)(double*,double*,double*,double*,double*);
 
34
extern void F_FUNC(e1xb,E1XB)(double*,double*);
 
35
extern void F_FUNC(e1z,E1Z)(Py_complex*,Py_complex*);
 
36
extern void F_FUNC(eix,EIX)(double*,double*);
 
37
extern void F_FUNC(cerror,CERROR)(Py_complex*,Py_complex*);
 
38
extern void F_FUNC(stvh0,STVH0)(double*,double*);
 
39
extern void F_FUNC(stvh1,STVH1)(double*,double*);
 
40
extern void F_FUNC(stvhv,STVHV)(double*,double*,double*);
 
41
extern void F_FUNC(stvl0,STVL0)(double*,double*);
 
42
extern void F_FUNC(stvl1,STVL1)(double*,double*);
 
43
extern void F_FUNC(stvlv,STVLV)(double*,double*,double*);
 
44
extern void F_FUNC(itsh0,ITSH0)(double*,double*);
 
45
extern void F_FUNC(itth0,ITTH0)(double*,double*);
 
46
extern void F_FUNC(itsl0,ITSL0)(double*,double*);
 
47
extern void F_FUNC(klvna,KLVNA)(double*,double*,double*,double*,double*,double*,double*,double*,double*);
 
48
extern void F_FUNC(itjya,ITJYA)(double*,double*,double*);
 
49
extern void F_FUNC(ittjya,ITTJYA)(double*,double*,double*);
 
50
extern void F_FUNC(itika,ITIKA)(double*,double*,double*);
 
51
extern void F_FUNC(ittika,ITTIKA)(double*,double*,double*);
 
52
extern void F_FUNC(cfc,CFC)(Py_complex*,Py_complex*,Py_complex*);
 
53
extern void F_FUNC(cfs,CFS)(Py_complex*,Py_complex*,Py_complex*);
 
54
extern void F_FUNC(cva2,CVA2)(int*,int*,double*,double*);
 
55
extern void F_FUNC(mtu0,MTU0)(int*,int*,double*,double*,double*,double*);
 
56
extern void F_FUNC(mtu12,MTU12)(int*,int*,int*,double*,double*,double*,double*,double*,double*);
 
57
extern void F_FUNC(lpmv,LPMV)(double*,int*,double*,double*);
 
58
extern void F_FUNC(pbwa,PBWA)(double*,double*,double*,double*,double*,double*);
 
59
extern void F_FUNC(pbdv,PBDV)(double*,double*,double*,double*,double*,double*);
 
60
extern void F_FUNC(pbvv,PBVV)(double*,double*,double*,double*,double*,double*);
 
61
extern void F_FUNC(segv,SEGV)(int*,int*,double*,int*,double*,double*);
 
62
extern void F_FUNC(aswfa,ASWFA)(int*,int*,double*,double*,int*,double*,double*,double*);
 
63
extern void F_FUNC(rswfp,RSWFP)(int*,int*,double*,double*,double*,int*,double*,double*,double*,double*);
 
64
extern void F_FUNC(rswfo,RSWFO)(int*,int*,double*,double*,double*,int*,double*,double*,double*,double*);
 
65
extern void F_FUNC(ffk,FFK)(int*,double*,double*,double*,double*,double*,double*,double*,double*,double*);
 
66
 
 
67
 
 
68
/* This must be linked with fortran
 
69
 */
 
70
 
 
71
Py_complex cgamma_wrap( Py_complex z) {
 
72
  int kf = 1;
 
73
  Py_complex cy;
 
74
 
 
75
  F_FUNC(cgama,CGAMA)(CADDR(z), &kf, CADDR(cy));
 
76
  return cy;
 
77
}
 
78
 
 
79
Py_complex clngamma_wrap( Py_complex z) {
 
80
  int kf = 0;
 
81
  Py_complex cy;
 
82
 
 
83
  F_FUNC(cgama,CGAMA)(CADDR(z), &kf, CADDR(cy));
 
84
  return cy;
 
85
}
 
86
 
 
87
Py_complex cpsi_wrap( Py_complex z) {
 
88
  Py_complex cy;
 
89
  
 
90
  if (IMAG(z)==0.0) {
 
91
    REAL(cy) = psi(REAL(z));
 
92
    IMAG(cy) = 0.0;
 
93
  }
 
94
  else {
 
95
    F_FUNC(cpsi,CPSI)(CADDR(z), CADDR(cy));
 
96
  }
 
97
  return cy;
 
98
}
 
99
 
 
100
Py_complex crgamma_wrap( Py_complex z) {
 
101
  int kf = 1;
 
102
  Py_complex cy;
 
103
  Py_complex cy2;
 
104
  double magsq;
 
105
 
 
106
  F_FUNC(cgama,CGAMA)(CADDR(z), &kf, CADDR(cy));
 
107
  magsq = ABSQ(cy);
 
108
  REAL(cy2) = REAL(cy) / magsq;
 
109
  IMAG(cy2) = -IMAG(cy) / magsq;
 
110
  return cy2;
 
111
}
 
112
 
 
113
Py_complex chyp2f1_wrap( double a, double b, double c, Py_complex z) {
 
114
  Py_complex outz;
 
115
  int l1, l0;
 
116
 
 
117
 
 
118
  l0 = ((c == floor(c)) && (c < 0));
 
119
  l1 = ((fabs(1-REAL(z)) < 1e-15) && (IMAG(z) == 0) && (c-a-b <= 0));
 
120
  if (l0 || l1) {
 
121
    REAL(outz) = INFINITY;
 
122
    IMAG(outz) = 0.0;
 
123
    return outz;
 
124
  }
 
125
  F_FUNC(hygfz, HYGFZ)(&a, &b, &c, &z, &outz);
 
126
  return outz;
 
127
}
 
128
 
 
129
Py_complex chyp1f1_wrap(double a, double b, Py_complex z) {
 
130
  Py_complex outz;
 
131
 
 
132
  F_FUNC(cchg,CCHG)(&a, &b, &z, &outz);
 
133
  if (REAL(outz) == 1e300) {
 
134
    REAL(outz) = INFINITY;
 
135
  }
 
136
  return outz;
 
137
}
 
138
 
 
139
 
 
140
double hypU_wrap(double a, double b, double x) {
 
141
  double out;
 
142
  int md; /* method code --- not returned */
 
143
 
 
144
  F_FUNC(chgu,CHGU)(&a, &b, &x, &out, &md);
 
145
  if (out == 1e300) out = INFINITY;
 
146
  return out;
 
147
  
 
148
}
 
149
 
 
150
 
 
151
int itairy_wrap(double x, double *apt, double *bpt, double *ant, double *bnt) {
 
152
  double tmp; 
 
153
  int flag = 0;
 
154
  
 
155
  if (x < 0) {
 
156
    x = -x;
 
157
    flag = 1;
 
158
  }
 
159
  F_FUNC(itairy,ITAIRY)(&x, apt, bpt, ant, bnt);
 
160
  if (flag) {  /* negative limit -- switch signs and roles */
 
161
    tmp = *apt;
 
162
    *apt = -*ant;
 
163
    *ant = -tmp;
 
164
    tmp = *bpt;
 
165
    *bpt = -*bnt;
 
166
    *bnt = -tmp;
 
167
  }
 
168
  return 0;
 
169
}
 
170
 
 
171
 
 
172
double exp1_wrap(double x) {
 
173
  double out;
 
174
  
 
175
  F_FUNC(e1xb,E1XB)(&x, &out);
 
176
  CONVINF(out);
 
177
  return out;
 
178
}
 
179
 
 
180
Py_complex cexp1_wrap(Py_complex z) {
 
181
  Py_complex outz;
 
182
  
 
183
  F_FUNC(e1z,E1Z)(&z, &outz);
 
184
  ZCONVINF(outz);
 
185
  return outz;
 
186
}
 
187
 
 
188
double expi_wrap(double x) {
 
189
  double out;
 
190
  
 
191
  F_FUNC(eix,EIX)(&x, &out);
 
192
  CONVINF(out);
 
193
  return out;
 
194
}
 
195
 
 
196
Py_complex cerf_wrap(Py_complex z) {
 
197
  Py_complex outz;
 
198
  
 
199
  F_FUNC(cerror,CERROR)(&z, &outz);
 
200
  return outz;
 
201
}
 
202
 
 
203
double struve_wrap(double v, double x) {
 
204
  double out;
 
205
  int flag=0;
 
206
 
 
207
  if ((v<-8.0) || (v>12.5)) {
 
208
    return struve(v, x);  /* from cephes */
 
209
  }
 
210
  if (v==0.0) {
 
211
    if (x < 0) {x = -x; flag=1;}
 
212
    F_FUNC(stvh0,STVH0)(&x,&out);
 
213
    CONVINF(out);
 
214
    if (flag) out = -out;
 
215
    return out;
 
216
  }
 
217
  if (v==1.0) {
 
218
    if (x < 0) x=-x;
 
219
    F_FUNC(stvh1,STVH1)(&x,&out);
 
220
    CONVINF(out);
 
221
    return out;
 
222
  }
 
223
  F_FUNC(stvhv,STVHV)(&v,&x,&out);
 
224
  CONVINF(out);
 
225
  return out;  
 
226
}
 
227
 
 
228
double modstruve_wrap(double v, double x) {
 
229
  double out;
 
230
  int flag=0;
 
231
 
 
232
  if ((x < 0) & (floor(v)!=v)) return NAN;
 
233
  if (v==0.0) {
 
234
    if (x < 0) {x = -x; flag=1;}
 
235
    F_FUNC(stvl0,STVl0)(&x,&out);
 
236
    CONVINF(out);
 
237
    if (flag) out = -out;
 
238
    return out;
 
239
  }
 
240
  if (v==1.0) {
 
241
    if (x < 0) x=-x;
 
242
    F_FUNC(stvl1,STVl1)(&x,&out);
 
243
    CONVINF(out);
 
244
    return out;
 
245
  }
 
246
  if (x<0) {
 
247
    x = -x;
 
248
    flag = 1;
 
249
  }
 
250
  F_FUNC(stvlv,STVlV)(&v,&x,&out);
 
251
  CONVINF(out);
 
252
  if (flag && (!((int)floor(v) % 2))) out = -out;
 
253
  return out;  
 
254
}
 
255
 
 
256
double itstruve0_wrap(double x) {
 
257
  double out;
 
258
 
 
259
  if (x<0) x=-x;
 
260
  F_FUNC(itsh0,ITSH0)(&x,&out);
 
261
  CONVINF(out);
 
262
  return out;
 
263
}
 
264
 
 
265
double it2struve0_wrap(double x) {
 
266
  double out;
 
267
  int flag=0;
 
268
  
 
269
  if (x<0) {x=-x; flag=1;}
 
270
  F_FUNC(itth0,ITTH0)(&x,&out);
 
271
  CONVINF(out);
 
272
  if (flag) {
 
273
    out = PI - out;
 
274
  }
 
275
  return out;
 
276
}
 
277
 
 
278
double itmodstruve0_wrap(double x) {
 
279
  double out;
 
280
 
 
281
  if (x<0) x=-x;
 
282
  F_FUNC(itsl0,ITSL0)(&x,&out);
 
283
  CONVINF(out);
 
284
  return out;
 
285
}
 
286
 
 
287
 
 
288
double ber_wrap(double x)
 
289
{
 
290
  Py_complex Be, Ke, Bep, Kep;
 
291
 
 
292
  if (x<0) x=-x;
 
293
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
294
  ZCONVINF(Be);
 
295
  return REAL(Be);
 
296
}
 
297
 
 
298
double bei_wrap(double x)
 
299
{
 
300
  Py_complex Be, Ke, Bep, Kep;
 
301
 
 
302
  if (x<0) x=-x;
 
303
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
304
  ZCONVINF(Be);
 
305
  return IMAG(Be);
 
306
}
 
307
 
 
308
double ker_wrap(double x)
 
309
{
 
310
  Py_complex Be, Ke, Bep, Kep;
 
311
 
 
312
  if (x<0) return NAN;
 
313
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
314
  ZCONVINF(Ke);
 
315
  return REAL(Ke);  
 
316
}
 
317
 
 
318
double kei_wrap(double x)
 
319
{
 
320
  Py_complex Be, Ke, Bep, Kep;
 
321
 
 
322
  if (x<0) return NAN;
 
323
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
324
  ZCONVINF(Ke);
 
325
  return IMAG(Ke);  
 
326
}
 
327
 
 
328
double berp_wrap(double x)
 
329
{
 
330
  Py_complex Be, Ke, Bep, Kep;
 
331
  int flag = 0;
 
332
 
 
333
  if (x<0) {x=-x; flag=1;}
 
334
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
335
  ZCONVINF(Bep);
 
336
  if (flag) return -REAL(Bep);
 
337
  return REAL(Bep);
 
338
}
 
339
 
 
340
double beip_wrap(double x)
 
341
{
 
342
  Py_complex Be, Ke, Bep, Kep;
 
343
  int flag = 0;
 
344
 
 
345
  if (x<0) {x=-x; flag=1;}
 
346
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
347
  ZCONVINF(Bep);
 
348
  if (flag) return -IMAG(Bep);
 
349
  return IMAG(Bep);
 
350
}
 
351
 
 
352
double kerp_wrap(double x)
 
353
{
 
354
  Py_complex Be, Ke, Bep, Kep;
 
355
 
 
356
  if (x<0) return NAN;
 
357
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
358
  ZCONVINF(Kep);
 
359
  return REAL(Kep);  
 
360
}
 
361
 
 
362
double keip_wrap(double x)
 
363
{
 
364
  Py_complex Be, Ke, Bep, Kep;
 
365
 
 
366
  if (x<0) return NAN;
 
367
  F_FUNC(klvna,KLVNA)(&x, CADDR(Be), CADDR(Ke), CADDR(Bep), CADDR(Kep));
 
368
  ZCONVINF(Kep);
 
369
  return IMAG(Kep);  
 
370
}
 
371
 
 
372
 
 
373
int kelvin_wrap(double x, Py_complex *Be, Py_complex *Ke, Py_complex *Bep, Py_complex *Kep) {
 
374
  int flag = 0;
 
375
  
 
376
  if (x<0) {x=-x; flag=1;}
 
377
  F_FUNC(klvna,KLVNA)(&x, F2C_CST(Be), F2C_CST(Ke), F2C_CST(Bep), F2C_CST(Kep));
 
378
  ZCONVINF(*Be);
 
379
  ZCONVINF(*Ke);
 
380
  ZCONVINF(*Bep);
 
381
  ZCONVINF(*Kep);
 
382
  if (flag) {
 
383
    REAL(*Bep) = -REAL(*Bep);
 
384
    IMAG(*Bep) = -IMAG(*Bep);
 
385
    REAL(*Ke) = NAN;
 
386
    IMAG(*Ke) = NAN;
 
387
    REAL(*Kep) = NAN;
 
388
    IMAG(*Kep) = NAN;    
 
389
  }    
 
390
  return 0;
 
391
}
 
392
 
 
393
/* Integrals of bessel functions */
 
394
 
 
395
/* int(j0(t),t=0..x) */
 
396
/* int(y0(t),t=0..x) */
 
397
 
 
398
int it1j0y0_wrap(double x, double *j0int, double *y0int)
 
399
{
 
400
  int flag = 0;
 
401
 
 
402
  if (x < 0) {x = -x; flag=1;}
 
403
  F_FUNC(itjya, ITJYA)(&x, j0int, y0int);
 
404
  if (flag) {
 
405
    *j0int = -(*j0int);
 
406
    *y0int = NAN;    /* domain error */
 
407
  }
 
408
  return 0;
 
409
}
 
410
 
 
411
/* int((1-j0(t))/t,t=0..x) */
 
412
/* int(y0(t)/t,t=x..inf) */
 
413
 
 
414
int it2j0y0_wrap(double x, double *j0int, double *y0int) 
 
415
{
 
416
  int flag = 0;
 
417
 
 
418
  if (x < 0) {x=-x; flag=1;}
 
419
  F_FUNC(ittjya, ITTJYA)(&x, j0int, y0int);
 
420
  if (flag) {
 
421
    *y0int = NAN;  /* domain error */
 
422
  }
 
423
  return 0;
 
424
}
 
425
 
 
426
/* Integrals of modified bessel functions */
 
427
 
 
428
int it1i0k0_wrap(double x, double *i0int, double *k0int)
 
429
{
 
430
  int flag = 0;
 
431
 
 
432
  if (x < 0) {x = -x; flag=1;}
 
433
  F_FUNC(itika, ITIKA)(&x, i0int, k0int);
 
434
  if (flag) {
 
435
    *i0int = -(*i0int);
 
436
    *k0int = NAN;    /* domain error */
 
437
  }
 
438
  return 0;
 
439
}
 
440
 
 
441
int it2i0k0_wrap(double x, double *i0int, double *k0int) 
 
442
{
 
443
  int flag = 0;
 
444
 
 
445
  if (x < 0) {x=-x; flag=1;}
 
446
  F_FUNC(ittika, ITTIKA)(&x, i0int, k0int);
 
447
  if (flag) {
 
448
    *k0int = NAN;  /* domain error */
 
449
  }
 
450
  return 0;
 
451
}
 
452
 
 
453
 
 
454
/* Fresnel integrals of complex numbers */
 
455
 
 
456
int cfresnl_wrap(Py_complex z, Py_complex *zfs, Py_complex *zfc)
 
457
{
 
458
  Py_complex zfd;
 
459
  F_FUNC(cfs,CFS)(&z,zfs,&zfd);
 
460
  F_FUNC(cfc,CFC)(&z,zfc,&zfd); 
 
461
  return 0;
 
462
}
 
463
 
 
464
/* Mathieu functions */
 
465
/* Characteristic values */
 
466
double cem_cva_wrap(double m, double q) {
 
467
  int int_m, kd=1;
 
468
  double out;
 
469
 
 
470
  if ((m < 0) || (m != floor(m))) 
 
471
    return NAN;
 
472
  int_m = (int )m;
 
473
  if (int_m % 2) kd=2;
 
474
  F_FUNC(cva2,CVA2)(&kd, &int_m, &q, &out);
 
475
  return out;               
 
476
}
 
477
 
 
478
double sem_cva_wrap(double m, double q) {
 
479
  int int_m, kd=4;
 
480
  double out;
 
481
 
 
482
  if ((m < 1) || (m != floor(m))) 
 
483
    return NAN;
 
484
  int_m = (int )m;
 
485
  if (int_m % 2) kd=3;
 
486
  F_FUNC(cva2,CVA2)(&kd, &int_m, &q, &out);
 
487
  return out;               
 
488
}
 
489
 
 
490
/* Mathieu functions */
 
491
int cem_wrap(double m, double q, double x, double *csf, double *csd)
 
492
{
 
493
  int int_m, kf=1;
 
494
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
495
    *csf = NAN;
 
496
    *csd = NAN;
 
497
  }
 
498
  int_m = (int )m;
 
499
  F_FUNC(mtu0,MTU0)(&kf,&int_m, &q, &x, csf, csd);
 
500
  return 0;  
 
501
}
 
502
 
 
503
int sem_wrap(double m, double q, double x, double *csf, double *csd)
 
504
{
 
505
  int int_m, kf=2;
 
506
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
507
    *csf = NAN;
 
508
    *csd = NAN;
 
509
  }
 
510
  int_m = (int )m;
 
511
  F_FUNC(mtu0,MTU0)(&kf,&int_m, &q, &x, csf, csd);
 
512
  return 0;  
 
513
}
 
514
 
 
515
 
 
516
int mcm1_wrap(double m, double q, double x, double *f1r, double *d1r)
 
517
{
 
518
  int int_m, kf=1, kc=1;
 
519
  double f2r, d2r;
 
520
 
 
521
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
522
    *f1r = NAN;
 
523
    *d1r = NAN;
 
524
  }
 
525
  int_m = (int )m;
 
526
  F_FUNC(mtu12,MTU12)(&kf,&kc,&int_m, &q, &x, f1r, d1r, &f2r, &d2r);
 
527
  return 0;  
 
528
}
 
529
 
 
530
int msm1_wrap(double m, double q, double x, double *f1r, double *d1r)
 
531
{
 
532
  int int_m, kf=2, kc=1;
 
533
  double f2r, d2r;
 
534
 
 
535
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
536
    *f1r = NAN;
 
537
    *d1r = NAN;
 
538
  }
 
539
  int_m = (int )m;
 
540
  F_FUNC(mtu12,MTU12)(&kf,&kc,&int_m, &q, &x, f1r, d1r, &f2r, &d2r);
 
541
  return 0;  
 
542
}
 
543
 
 
544
int mcm2_wrap(double m, double q, double x, double *f2r, double *d2r)
 
545
{
 
546
  int int_m, kf=1, kc=2;
 
547
  double f1r, d1r;
 
548
 
 
549
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
550
    *f2r = NAN;
 
551
    *d2r = NAN;
 
552
  }
 
553
  int_m = (int )m;
 
554
  F_FUNC(mtu12,MTU12)(&kf,&kc,&int_m, &q, &x, &f1r, &d1r, f2r, d2r);
 
555
  return 0;  
 
556
}
 
557
 
 
558
int msm2_wrap(double m, double q, double x, double *f2r, double *d2r)
 
559
{
 
560
  int int_m, kf=2, kc=2;
 
561
  double f1r, d1r;
 
562
 
 
563
  if ((m < 1) || (m != floor(m)) || (q<0)) {
 
564
    *f2r = NAN;
 
565
    *d2r = NAN;
 
566
  }
 
567
  int_m = (int )m;
 
568
  F_FUNC(mtu12,MTU12)(&kf,&kc,&int_m, &q, &x, &f1r, &d1r, f2r, d2r);
 
569
  return 0;  
 
570
}
 
571
 
 
572
 
 
573
double pmv_wrap(double m, double v, double x){
 
574
  int int_m;
 
575
  double out;
 
576
 
 
577
  if (m != floor(m)) return NAN;
 
578
  int_m = (int ) m;
 
579
  F_FUNC(lpmv,LPMV)(&v, &int_m, &x, &out);
 
580
  return out;
 
581
}
 
582
 
 
583
 
 
584
/* if x > 0 return w1f and w1d.
 
585
    otherwise return w2f and w2d (after abs(x))
 
586
*/
 
587
int pbwa_wrap(double a, double x, double *wf, double *wd) {
 
588
  int flag = 0;
 
589
  double w1f, w1d, w2f, w2d;
 
590
   
 
591
  if (x < 0) {x=-x; flag=1;}
 
592
  F_FUNC(pbwa,PBWA)(&a, &x, &w1f, &w1d, &w2f, &w2d);
 
593
  if (flag) {
 
594
    *wf = w2f;
 
595
    *wd = w2d;
 
596
  }
 
597
  else {
 
598
    *wf = w1f;
 
599
    *wd = w1d;
 
600
  }
 
601
  return 0;
 
602
}
 
603
 
 
604
int pbdv_wrap(double v, double x, double *pdf, double *pdd) {
 
605
 
 
606
  double *dv;
 
607
  double *dp;
 
608
  int num;
 
609
 
 
610
  num = ABS((int) v)+1;    
 
611
  dv = (double *)PyMem_Malloc(sizeof(double)*2*num);
 
612
  if (dv==NULL) {
 
613
    printf("Warning: Memory allocation error.\n"); 
 
614
    *pdf = NAN;
 
615
    *pdd = NAN;
 
616
    return -1;
 
617
  }
 
618
  dp = dv + num;
 
619
  F_FUNC(pbdv,PBDV)(&v, &x, dv, dp, pdf, pdd);
 
620
  PyMem_Free(dv);
 
621
  return 0;
 
622
}
 
623
 
 
624
int pbvv_wrap(double v, double x, double *pvf, double *pvd) {
 
625
  double *vv;
 
626
  double *vp;
 
627
  int num;
 
628
 
 
629
  num = ABS((int) v)+1;    
 
630
  vv = (double *)PyMem_Malloc(sizeof(double)*2*num);
 
631
  if (vv==NULL) {
 
632
    printf("Warning: Memory allocation error.\n"); 
 
633
    *pvf = NAN;
 
634
    *pvd = NAN;
 
635
    return -1;
 
636
  }
 
637
  vp = vv + num;
 
638
  F_FUNC(pbvv,PBVV)(&v, &x, vv, vp, pvf, pvd);
 
639
  PyMem_Free(vv);
 
640
  return 0;
 
641
}
 
642
 
 
643
double prolate_segv_wrap(double m, double n, double c)
 
644
{
 
645
  int kd=1;
 
646
  int int_m, int_n;
 
647
  double cv, *eg;
 
648
 
 
649
  if ((m<0) || (n<m) || (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
650
    return NAN;
 
651
  }
 
652
  int_m = (int) m;
 
653
  int_n = (int) n;
 
654
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
655
  if (eg==NULL) {
 
656
    printf("Warning: Memory allocation error.\n"); 
 
657
    return NAN;
 
658
  }
 
659
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
660
  PyMem_Free(eg);
 
661
  return cv;
 
662
}
 
663
 
 
664
double oblate_segv_wrap(double m, double n, double c)
 
665
{
 
666
  int kd=-1;
 
667
  int int_m, int_n;
 
668
  double cv, *eg;
 
669
 
 
670
  if ((m<0) || (n<m) || (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
671
    return NAN;
 
672
  }
 
673
  int_m = (int) m;
 
674
  int_n = (int) n;
 
675
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
676
  if (eg==NULL) {
 
677
    printf("Warning: Memory allocation error.\n"); 
 
678
    return NAN;
 
679
  }
 
680
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
681
  PyMem_Free(eg);
 
682
  return cv;
 
683
}
 
684
 
 
685
 
 
686
double prolate_aswfa_nocv_wrap(double m, double n, double c, double x, double *s1d)
 
687
{
 
688
  int kd = 1;
 
689
  int int_m, int_n;
 
690
  double cv, s1f, *eg;
 
691
 
 
692
  if ((x >=1) || (x <=-1) || (m<0) || (n<m) || \
 
693
      (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
694
    *s1d = NAN;
 
695
    return NAN;
 
696
  }
 
697
  int_m = (int )m;
 
698
  int_n = (int )n;
 
699
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
700
  if (eg==NULL) {
 
701
    printf("Warning: Memory allocation error.\n"); 
 
702
    *s1d = NAN;
 
703
    return NAN;
 
704
  }
 
705
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
706
  F_FUNC(aswfa,ASWFA)(&int_m,&int_n,&c,&x,&kd,&cv,&s1f,s1d);
 
707
  PyMem_Free(eg);
 
708
  return s1f;
 
709
}
 
710
 
 
711
 
 
712
double oblate_aswfa_nocv_wrap(double m, double n, double c, double x, double *s1d)
 
713
{
 
714
  int kd = -1;
 
715
  int int_m, int_n;
 
716
  double cv, s1f, *eg;
 
717
 
 
718
  if ((x >=1) || (x <=-1) || (m<0) || (n<m) || \
 
719
      (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
720
    *s1d = NAN;
 
721
    return NAN;
 
722
  }
 
723
  int_m = (int )m;
 
724
  int_n = (int )n;
 
725
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
726
  if (eg==NULL) {
 
727
    printf("Warning: Memory allocation error.\n"); 
 
728
    *s1d = NAN;
 
729
    return NAN;
 
730
  }
 
731
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
732
  F_FUNC(aswfa,ASWFA)(&int_m,&int_n,&c,&x,&kd,&cv,&s1f,s1d);
 
733
  PyMem_Free(eg);
 
734
  return s1f;
 
735
}
 
736
 
 
737
 
 
738
int prolate_aswfa_wrap(double m, double n, double c, double cv, double x, double *s1f, double *s1d)
 
739
{
 
740
  int kd = 1;
 
741
  int int_m, int_n;
 
742
 
 
743
  if ((x >=1) || (x <=-1) || (m<0) || (n<m) || \
 
744
      (m!=floor(m)) || (n!=floor(n))) {
 
745
    *s1f = NAN;
 
746
    *s1d = NAN;
 
747
    return 0;
 
748
  }
 
749
  int_m = (int )m;
 
750
  int_n = (int )n;
 
751
  F_FUNC(aswfa,ASWFA)(&int_m,&int_n,&c,&x,&kd,&cv,s1f,s1d);
 
752
  return 0;
 
753
}
 
754
 
 
755
 
 
756
int oblate_aswfa_wrap(double m, double n, double c, double cv, double x, double *s1f, double *s1d)
 
757
{
 
758
  int kd = -1;
 
759
  int int_m, int_n;
 
760
 
 
761
  if ((x >=1) || (x <=-1) || (m<0) || (n<m) || \
 
762
      (m!=floor(m)) || (n!=floor(n))) {
 
763
    *s1f = NAN;
 
764
    *s1d = NAN;
 
765
    return 0;
 
766
  }
 
767
  int_m = (int )m;
 
768
  int_n = (int )n;
 
769
  F_FUNC(aswfa,ASWFA)(&int_m,&int_n,&c,&x,&kd,&cv,s1f,s1d);
 
770
  return 0;
 
771
}
 
772
 
 
773
 
 
774
double prolate_radial1_nocv_wrap(double m, double n, double c, double x, double *r1d)
 
775
{
 
776
  int kf=1, kd=1;
 
777
  double r2f, r2d, r1f, cv, *eg;
 
778
  int int_m, int_n;
 
779
 
 
780
  if ((x <=1.0) || (m<0) || (n<m) || \
 
781
     (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
782
    *r1d = NAN;
 
783
    return NAN;
 
784
  }
 
785
  int_m = (int )m;
 
786
  int_n = (int )n;
 
787
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
788
  if (eg==NULL) {
 
789
    printf("Warning: Memory allocation error.\n"); 
 
790
    *r1d = NAN;
 
791
    return NAN;
 
792
  }
 
793
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
794
  F_FUNC(rswfp,RSWFP)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,r1d,&r2f,&r2d);
 
795
  PyMem_Free(eg);
 
796
  return r1f;
 
797
}
 
798
 
 
799
double prolate_radial2_nocv_wrap(double m, double n, double c, double x, double *r2d)
 
800
{
 
801
  int kf=2, kd=1;
 
802
  double r1f, r1d, r2f, cv, *eg;
 
803
  int int_m, int_n;
 
804
 
 
805
  if ((x <=1.0) || (m<0) || (n<m) || \
 
806
     (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
807
    *r2d = NAN;
 
808
    return NAN;
 
809
  }
 
810
  int_m = (int )m;
 
811
  int_n = (int )n;
 
812
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
813
  if (eg==NULL) {
 
814
    printf("Warning: Memory allocation error.\n"); 
 
815
    *r2d = NAN;
 
816
    return NAN;
 
817
  }
 
818
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
819
  F_FUNC(rswfp,RSWFP)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,&r1d,&r2f,r2d);
 
820
  PyMem_Free(eg);
 
821
  return r2f;
 
822
}
 
823
 
 
824
int prolate_radial1_wrap(double m, double n, double c, double cv, double x, double *r1f, double *r1d)
 
825
{
 
826
  int kf=1;
 
827
  double r2f, r2d;
 
828
  int int_m, int_n;
 
829
 
 
830
  if ((x <= 1.0) || (m<0) || (n<m) || \
 
831
     (m!=floor(m)) || (n!=floor(n))) {
 
832
    *r1f = NAN;
 
833
    *r1d = NAN;
 
834
    return 0;
 
835
  }
 
836
  int_m = (int )m;
 
837
  int_n = (int )n;
 
838
  F_FUNC(rswfp,RSWFP)(&int_m,&int_n,&c,&x,&cv,&kf,r1f,r1d,&r2f,&r2d);
 
839
  return 0;  
 
840
}
 
841
 
 
842
int prolate_radial2_wrap(double m, double n, double c, double cv, double x, double *r2f, double *r2d)
 
843
{
 
844
  int kf=2;
 
845
  double r1f, r1d;
 
846
  int int_m, int_n;
 
847
 
 
848
  if ((x <= 1.0) || (m<0) || (n<m) || \
 
849
     (m!=floor(m)) || (n!=floor(n))) {
 
850
    *r2f = NAN;
 
851
    *r2d = NAN;
 
852
    return 0;
 
853
  }
 
854
  int_m = (int )m;
 
855
  int_n = (int )n;
 
856
  F_FUNC(rswfp,RSWFP)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,&r1d,r2f,r2d);
 
857
  return 0;  
 
858
}
 
859
 
 
860
double oblate_radial1_nocv_wrap(double m, double n, double c, double x, double *r1d)
 
861
{
 
862
  int kf=1, kd=-1;
 
863
  double r2f, r2d, r1f, cv, *eg;
 
864
  int int_m, int_n;
 
865
 
 
866
  if ((x < 0.0) || (m<0) || (n<m) || \
 
867
     (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
868
    *r1d = NAN;
 
869
    return NAN;
 
870
  }
 
871
  int_m = (int )m;
 
872
  int_n = (int )n;
 
873
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
874
  if (eg==NULL) {
 
875
    printf("Warning: Memory allocation error.\n"); 
 
876
    *r1d = NAN;
 
877
    return NAN;
 
878
  }
 
879
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
880
  F_FUNC(rswfo,RSWFO)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,r1d,&r2f,&r2d);
 
881
  PyMem_Free(eg);
 
882
  return r1f;
 
883
}
 
884
 
 
885
double oblate_radial2_nocv_wrap(double m, double n, double c, double x, double *r2d)
 
886
{
 
887
  int kf=2, kd=-1;
 
888
  double r1f, r1d, r2f, cv, *eg;
 
889
  int int_m, int_n;
 
890
 
 
891
  if ((x < 0.0) || (m<0) || (n<m) || \
 
892
     (m!=floor(m)) || (n!=floor(n)) || ((n-m)>198)) {
 
893
    *r2d = NAN;
 
894
    return NAN;
 
895
  }
 
896
  int_m = (int )m;
 
897
  int_n = (int )n;
 
898
  eg = (double *)PyMem_Malloc(sizeof(double)*(n-m+2));
 
899
  if (eg==NULL) {
 
900
    printf("Warning: Memory allocation error.\n"); 
 
901
    *r2d = NAN;
 
902
    return NAN;
 
903
  }
 
904
  F_FUNC(segv,SEGV)(&int_m,&int_n,&c,&kd,&cv,eg);
 
905
  F_FUNC(rswfo,RSWFO)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,&r1d,&r2f,r2d);
 
906
  PyMem_Free(eg);
 
907
  return r2f;
 
908
}
 
909
 
 
910
int oblate_radial1_wrap(double m, double n, double c, double cv, double x, double *r1f, double *r1d)
 
911
{
 
912
  int kf=1;
 
913
  double r2f, r2d;
 
914
  int int_m, int_n;
 
915
 
 
916
  if ((x <0.0) || (m<0) || (n<m) || \
 
917
     (m!=floor(m)) || (n!=floor(n))) {
 
918
    *r1f = NAN;
 
919
    *r1d = NAN;
 
920
    return 0;
 
921
  }
 
922
  int_m = (int )m;
 
923
  int_n = (int )n;
 
924
  F_FUNC(rswfo,RSWFO)(&int_m,&int_n,&c,&x,&cv,&kf,r1f,r1d,&r2f,&r2d);
 
925
  return 0;  
 
926
}
 
927
 
 
928
int oblate_radial2_wrap(double m, double n, double c, double cv, double x, double *r2f, double *r2d)
 
929
{
 
930
  int kf=2;
 
931
  double r1f, r1d;
 
932
  int int_m, int_n;
 
933
 
 
934
  if ((x <0.0) || (m<0) || (n<m) || \
 
935
     (m!=floor(m)) || (n!=floor(n))) {
 
936
    *r2f = NAN;
 
937
    *r2d = NAN;
 
938
    return 0;
 
939
  }
 
940
  int_m = (int )m;
 
941
  int_n = (int )n;
 
942
  F_FUNC(rswfo,RSWFO)(&int_m,&int_n,&c,&x,&cv,&kf,&r1f,&r1d,r2f,r2d);
 
943
  return 0;  
 
944
}
 
945
 
 
946
 
 
947
int modified_fresnel_plus_wrap(double x, Py_complex *Fplus, Py_complex *Kplus)
 
948
{
 
949
  int ks=0;
 
950
  double fm, fa, gm, ga;
 
951
  
 
952
  F_FUNC(ffk,FFK)(&ks,&x,F2C_CST(Fplus),&fm,&fa,F2C_CST(Kplus),&gm,&ga);
 
953
  return 0;
 
954
}
 
955
 
 
956
int modified_fresnel_minus_wrap(double x, Py_complex *Fminus, Py_complex *Kminus)
 
957
{
 
958
  int ks=1;
 
959
  double fm, fa, gm, ga;
 
960
  
 
961
  F_FUNC(ffk,FFK)(&ks,&x,F2C_CST(Fminus),&fm,&fa,F2C_CST(Kminus),&gm,&ga);
 
962
  return 0;
 
963
}
 
964