libdspl-2.0
Digital Signal Processing Algorithm Library
ellipj.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of DSPL.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include "dspl.h"
27 
28 
29 
30 #ifdef DOXYGEN_ENGLISH
31 
64 #endif
65 #ifdef DOXYGEN_RUSSIAN
66 
102 #endif
103 int DSPL_API ellip_acd(double* w, int n, double k, double* u)
104 {
105  double lnd[ELLIP_ITER], t;
106  int i, m;
107 
108  if(!u || !w)
109  return ERROR_PTR;
110  if(n<1)
111  return ERROR_SIZE;
112  if(k < 0.0 || k>= 1.0)
113  return ERROR_ELLIP_MODULE;
114 
115  ellip_landen(k,ELLIP_ITER, lnd);
116 
117  for(m = 0; m < n; m++)
118  {
119  u[m] = w[m];
120  for(i = 1; i < ELLIP_ITER; i++)
121  {
122  t = lnd[i-1]*u[m];
123  t *= t;
124  t = 1.0 + sqrt(1.0 - t);
125  u[m] = 2.0 * u[m] / (t+t*lnd[i]);
126  }
127  u[m] = 2.0 * acos(u[m]) / M_PI;
128  }
129  return RES_OK;
130 }
131 
132 
133 
134 
135 
136 
137 #ifdef DOXYGEN_ENGLISH
138 
171 #endif
172 #ifdef DOXYGEN_RUSSIAN
173 
210 #endif
211 int DSPL_API ellip_acd_cmplx(complex_t* w, int n, double k, complex_t* u)
212 {
213  double lnd[ELLIP_ITER], t;
214  complex_t tmp0, tmp1;
215  int i, m;
216 
217  if(!u || !w)
218  return ERROR_PTR;
219  if(n<1)
220  return ERROR_SIZE;
221  if(k < 0.0 || k>= 1.0)
222  return ERROR_ELLIP_MODULE;
223 
224  ellip_landen(k,ELLIP_ITER, lnd);
225 
226  for(m = 0; m < n; m++)
227  {
228  RE(u[m]) = RE(w[m]);
229  IM(u[m]) = IM(w[m]);
230  for(i = 1; i < ELLIP_ITER; i++)
231  {
232  RE(tmp0) = lnd[i-1]*RE(u[m]);
233  IM(tmp0) = lnd[i-1]*IM(u[m]);
234  RE(tmp1) = 1.0 - CMRE(tmp0, tmp0);
235  IM(tmp1) = - CMIM(tmp0, tmp0);
236 
237  sqrt_cmplx(&tmp1, 1, &tmp0);
238  RE(tmp0) += 1.0;
239 
240  RE(tmp1) = RE(tmp0) * (1.0 + lnd[i]);
241  IM(tmp1) = IM(tmp0) * (1.0 + lnd[i]);
242 
243  t = 2.0 / ABSSQR(tmp1);
244 
245  RE(tmp0) = t * CMCONJRE(u[m], tmp1);
246  IM(tmp0) = t * CMCONJIM(u[m], tmp1);
247 
248  RE(u[m]) = RE(tmp0);
249  IM(u[m]) = IM(tmp0);
250  }
251  acos_cmplx(&tmp0, 1, u+m);
252  t = 2.0 / M_PI;
253  RE(u[m]) *= t;
254  IM(u[m]) *= t;
255  }
256  return RES_OK;
257 }
258 
259 
260 
261 
262 
263 #ifdef DOXYGEN_ENGLISH
264 
297 #endif
298 #ifdef DOXYGEN_RUSSIAN
299 
338 #endif
339 int DSPL_API ellip_asn(double* w, int n, double k, double* u)
340 {
341  double lnd[ELLIP_ITER], t;
342  int i, m;
343 
344  if(!u || !w)
345  return ERROR_PTR;
346  if(n<1)
347  return ERROR_SIZE;
348  if(k < 0.0 || k>= 1.0)
349  return ERROR_ELLIP_MODULE;
350 
351  ellip_landen(k,ELLIP_ITER, lnd);
352 
353  for(m = 0; m < n; m++)
354  {
355  u[m] = w[m];
356  for(i = 1; i < ELLIP_ITER; i++)
357  {
358  t = lnd[i-1]*u[m];
359  t *= t;
360  t = 1.0 + sqrt(1.0 - t);
361  u[m] = 2.0 * u[m] / (t+t*lnd[i]);
362  }
363  u[m] = 2.0 * asin(u[m]) / M_PI;
364  }
365  return RES_OK;
366 }
367 
368 
369 
370 
371 
372 #ifdef DOXYGEN_ENGLISH
373 
406 #endif
407 #ifdef DOXYGEN_RUSSIAN
408 
446 #endif
447 int DSPL_API ellip_asn_cmplx(complex_t* w, int n, double k, complex_t* u)
448 {
449  double lnd[ELLIP_ITER], t;
450  complex_t tmp0, tmp1;
451  int i, m;
452 
453  if(!u || !w)
454  return ERROR_PTR;
455  if(n<1)
456  return ERROR_SIZE;
457  if(k < 0.0 || k>= 1.0)
458  return ERROR_ELLIP_MODULE;
459 
460  ellip_landen(k,ELLIP_ITER, lnd);
461 
462  for(m = 0; m < n; m++)
463  {
464  RE(u[m]) = RE(w[m]);
465  IM(u[m]) = IM(w[m]);
466  for(i = 1; i < ELLIP_ITER; i++)
467  {
468  RE(tmp0) = lnd[i-1]*RE(u[m]);
469  IM(tmp0) = lnd[i-1]*IM(u[m]);
470  RE(tmp1) = 1.0 - CMRE(tmp0, tmp0);
471  IM(tmp1) = - CMIM(tmp0, tmp0);
472 
473  sqrt_cmplx(&tmp1, 1, &tmp0);
474  RE(tmp0) += 1.0;
475 
476  RE(tmp1) = RE(tmp0) * (1.0 + lnd[i]);
477  IM(tmp1) = IM(tmp0) * (1.0 + lnd[i]);
478 
479  t = 2.0 / ABSSQR(tmp1);
480 
481  RE(tmp0) = t * CMCONJRE(u[m], tmp1);
482  IM(tmp0) = t * CMCONJIM(u[m], tmp1);
483 
484  RE(u[m]) = RE(tmp0);
485  IM(u[m]) = IM(tmp0);
486  }
487 
488  asin_cmplx(&tmp0, 1, u+m);
489  t = 2.0 / M_PI;
490  RE(u[m]) *= t;
491  IM(u[m]) *= t;
492  }
493  return RES_OK;
494 }
495 
496 
497 
498 
499 #ifdef DOXYGEN_ENGLISH
500 
543 #endif
544 #ifdef DOXYGEN_RUSSIAN
545 
593 #endif
594 int DSPL_API ellip_cd(double* u, int n, double k, double* y)
595 {
596  double lnd[ELLIP_ITER];
597  int i, m;
598 
599  if(!u || !y)
600  return ERROR_PTR;
601  if(n<1)
602  return ERROR_SIZE;
603  if(k < 0.0 || k>= 1.0)
604  return ERROR_ELLIP_MODULE;
605 
606  ellip_landen(k,ELLIP_ITER, lnd);
607 
608  for(m = 0; m < n; m++)
609  {
610  y[m] = cos(u[m] * M_PI * 0.5);
611  for(i = ELLIP_ITER-1; i>0; i--)
612  {
613  y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]);
614  }
615  }
616  return RES_OK;
617 }
618 
619 
620 
621 
622 
623 
624 
625 #ifdef DOXYGEN_ENGLISH
626 
660 #endif
661 #ifdef DOXYGEN_RUSSIAN
662 
700 #endif
701 int DSPL_API ellip_cd_cmplx(complex_t* u, int n, double k, complex_t* y)
702 {
703  double lnd[ELLIP_ITER], t;
704  int i, m;
705  complex_t tmp;
706 
707  if(!u || !y)
708  return ERROR_PTR;
709  if(n<1)
710  return ERROR_SIZE;
711  if(k < 0.0 || k>= 1.0)
712  return ERROR_ELLIP_MODULE;
713 
714  ellip_landen(k,ELLIP_ITER, lnd);
715 
716  for(m = 0; m < n; m++)
717  {
718  RE(tmp) = RE(u[m]) * M_PI * 0.5;
719  IM(tmp) = IM(u[m]) * M_PI * 0.5;
720 
721  cos_cmplx(&tmp, 1, y+m);
722 
723  for(i = ELLIP_ITER-1; i>0; i--)
724  {
725  t = 1.0 / ABSSQR(y[m]);
726 
727  RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i];
728  IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i];
729 
730  t = (1.0 + lnd[i]) / ABSSQR(tmp);
731 
732  RE(y[m]) = RE(tmp) * t;
733  IM(y[m]) = -IM(tmp) * t;
734  }
735  }
736  return RES_OK;
737 }
738 
739 
740 
741 
742 
743 
744 #ifdef DOXYGEN_ENGLISH
745 
814 #endif
815 #ifdef DOXYGEN_RUSSIAN
816 
888 #endif
889 int DSPL_API ellip_landen(double k, int n, double* y)
890 {
891  int i;
892  y[0] = k;
893 
894  if(!y)
895  return ERROR_PTR;
896  if(n < 1)
897  return ERROR_SIZE;
898  if(k < 0.0 || k>= 1.0)
899  return ERROR_ELLIP_MODULE;
900 
901  for(i = 1; i < n; i++)
902  {
903  y[i] = y[i-1] / (1.0 + sqrt(1.0 - y[i-1] * y[i-1]));
904  y[i] *= y[i];
905  }
906 
907  return RES_OK;
908 }
909 
910 
911 
912 
913 
914 #ifdef DOXYGEN_ENGLISH
915 
916 #endif
917 #ifdef DOXYGEN_RUSSIAN
918 
919 #endif
920 int DSPL_API ellip_modulareq(double rp, double rs, int ord, double *k)
921 {
922  double ep, es, ke, kp, t, sn = 0.0;
923  int i, L, r;
924 
925  if(rp < 0 || rp == 0)
926  return ERROR_FILTER_RP;
927  if(rs < 0 || rs == 0)
928  return ERROR_FILTER_RS;
929  if(ord < 1)
930  return ERROR_FILTER_ORD;
931  if(!k)
932  return ERROR_PTR;
933 
934 
935  ep = sqrt(pow(10.0, rp*0.1)-1.0);
936  es = sqrt(pow(10.0, rs*0.1)-1.0);
937 
938  ke = ep/es;
939 
940  ke = sqrt(1.0 - ke*ke);
941 
942  r = ord % 2;
943  L = (ord-r)/2;
944 
945  kp = 1.0;
946  for(i = 0; i < L; i++)
947  {
948  t = (double)(2*i+1) / (double)ord;
949  ellip_sn(&t, 1, ke, &sn);
950  sn*=sn;
951  kp *= sn*sn;
952  }
953 
954  kp *= pow(ke, (double)ord);
955  *k = sqrt(1.0 - kp*kp);
956 
957  return RES_OK;
958 
959 }
960 
961 
962 
963 #ifdef DOXYGEN_ENGLISH
964 
965 #endif
966 #ifdef DOXYGEN_RUSSIAN
967 
968 #endif
969 int DSPL_API ellip_rat(double* w, int n, int ord, double k, double* u)
970 {
971  double t, xi, w2, xi2, k2;
972  int i, m, r, L;
973 
974  if(!u || !w)
975  return ERROR_PTR;
976  if(n<1)
977  return ERROR_SIZE;
978  if(k < 0.0 || k>= 1.0)
979  return ERROR_ELLIP_MODULE;
980 
981  r = ord%2;
982  L = (ord-r)/2;
983 
984  if(r)
985  memcpy(u, w, n*sizeof(double));
986  else
987  {
988  for(m = 0; m < n; m++)
989  {
990  u[m] = 1.0;
991  }
992  }
993 
994  k2 = k*k;
995  for(i = 0; i < L; i++)
996  {
997  t = (double)(2*i+1) / (double)ord;
998  ellip_cd(&t, 1, k, &xi);
999  xi2 = xi*xi;
1000  for(m = 0; m < n; m++)
1001  {
1002  w2 = w[m]*w[m];
1003  u[m] *= (w2 - xi2) / (1.0 - w2 * k2 * xi2);
1004  u[m] *= (1.0 - k2*xi2) / (1.0 - xi2);
1005  }
1006  }
1007  return RES_OK;
1008 }
1009 
1010 
1011 
1012 
1013 
1014 
1015 #ifdef DOXYGEN_ENGLISH
1016 
1059 #endif
1060 #ifdef DOXYGEN_RUSSIAN
1061 
1110 #endif
1111 int DSPL_API ellip_sn(double* u, int n, double k, double* y)
1112 {
1113  double lnd[ELLIP_ITER];
1114  int i, m;
1115 
1116  if(!u || !y)
1117  return ERROR_PTR;
1118  if(n<1)
1119  return ERROR_SIZE;
1120  if(k < 0.0 || k>= 1.0)
1121  return ERROR_ELLIP_MODULE;
1122 
1123  ellip_landen(k,ELLIP_ITER, lnd);
1124 
1125 
1126  for(m = 0; m < n; m++)
1127  {
1128  y[m] = sin(u[m] * M_PI * 0.5);
1129  for(i = ELLIP_ITER-1; i>0; i--)
1130  {
1131  y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]);
1132  }
1133  }
1134  return RES_OK;
1135 }
1136 
1137 
1138 
1139 
1140 #ifdef DOXYGEN_ENGLISH
1141 
1175 #endif
1176 #ifdef DOXYGEN_RUSSIAN
1177 
1215 #endif
1216 int DSPL_API ellip_sn_cmplx(complex_t* u, int n, double k, complex_t* y)
1217 {
1218  double lnd[ELLIP_ITER], t;
1219  int i, m;
1220  complex_t tmp;
1221 
1222  if(!u || !y)
1223  return ERROR_PTR;
1224  if(n<1)
1225  return ERROR_SIZE;
1226  if(k < 0.0 || k>= 1.0)
1227  return ERROR_ELLIP_MODULE;
1228 
1229  ellip_landen(k,ELLIP_ITER, lnd);
1230 
1231 
1232  for(m = 0; m < n; m++)
1233  {
1234  RE(tmp) = RE(u[m]) * M_PI * 0.5;
1235  IM(tmp) = IM(u[m]) * M_PI * 0.5;
1236 
1237  sin_cmplx(&tmp, 1, y+m);
1238 
1239  for(i = ELLIP_ITER-1; i>0; i--)
1240  {
1241  t = 1.0 / ABSSQR(y[m]);
1242 
1243  RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i];
1244  IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i];
1245 
1246  t = (1.0 + lnd[i]) / ABSSQR(tmp);
1247 
1248  RE(y[m]) = RE(tmp) * t;
1249  IM(y[m]) = -IM(tmp) * t;
1250 
1251  }
1252  }
1253  return RES_OK;
1254 }
1255 
int DSPL_API ellip_cd(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:594
#define ERROR_FILTER_ORD
The filter order is incorrect. The order of the filter must be specified by a positive integer value.
Definition: dspl.h:575
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:420
int DSPL_API acos_cmplx(complex_t *x, int n, complex_t *y)
The inverse of the cosine function the complex vector argument x.
Definition: math.c:141
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:610
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:618
#define ERROR_FILTER_RP
The filter non-uniformity parameter in the passband is set incorrectly. This parameter is specified i...
Definition: dspl.h:577
int DSPL_API ellip_acd(double *w, int n, double k, double *u)
Inverse Jacobi elliptic function of the real vector argument.
Definition: ellipj.c:103
#define ERROR_FILTER_RS
The filter suppression parameter in the boom bar is not set correctly. This parameter is specified in...
Definition: dspl.h:578
int DSPL_API ellip_sn_cmplx(complex_t *u, int n, double k, complex_t *y)
Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:1216
int DSPL_API ellip_sn(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:1111
int DSPL_API ellip_asn_cmplx(complex_t *w, int n, double k, complex_t *u)
Inverse Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:447
double complex_t[2]
Complex data type.
Definition: dspl.h:86
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:558
#define ABSSQR(x)
The macro returns the square of the modulus of a complex number x.
Definition: dspl.h:534
int DSPL_API sqrt_cmplx(complex_t *x, int n, complex_t *y)
Square root of the complex vector argguument x.
Definition: math.c:1307
int DSPL_API sin_cmplx(complex_t *x, int n, complex_t *y)
The sine function the complex vector argument x.
Definition: math.c:849
int DSPL_API ellip_asn(double *w, int n, double k, double *u)
Inverse Jacobi elliptic function of real vector argument.
Definition: ellipj.c:339
int DSPL_API ellip_cd_cmplx(complex_t *u, int n, double k, complex_t *y)
Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:701
#define ERROR_ELLIP_MODULE
The modulus of the Jacobi elliptic integral must be from 0 to 1. This error occurs when calculating J...
Definition: dspl.h:569
int DSPL_API asin_cmplx(complex_t *x, int n, complex_t *y)
The inverse of the sine function the complex vector argument x.
Definition: math.c:271
int DSPL_API cos_cmplx(complex_t *x, int n, complex_t *y)
The cosine function the complex vector argument x.
Definition: math.c:585
int DSPL_API ellip_acd_cmplx(complex_t *w, int n, double k, complex_t *u)
Inverse Jacobi elliptic function of complex vector argument.
Definition: ellipj.c:211
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:478
int DSPL_API ellip_landen(double k, int n, double *y)
Function calculates complete elliptical integral coefficients .
Definition: ellipj.c:889