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 
534 #endif
535 #ifdef DOXYGEN_RUSSIAN
536 
576 #endif
577 int DSPL_API ellip_cd(double* u, int n, double k, double* y)
578 {
579  double lnd[ELLIP_ITER];
580  int i, m;
581 
582  if(!u || !y)
583  return ERROR_PTR;
584  if(n<1)
585  return ERROR_SIZE;
586  if(k < 0.0 || k>= 1.0)
587  return ERROR_ELLIP_MODULE;
588 
589  ellip_landen(k,ELLIP_ITER, lnd);
590 
591  for(m = 0; m < n; m++)
592  {
593  y[m] = cos(u[m] * M_PI * 0.5);
594  for(i = ELLIP_ITER-1; i>0; i--)
595  {
596  y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]);
597  }
598  }
599  return RES_OK;
600 }
601 
602 
603 
604 
605 
606 
607 
608 #ifdef DOXYGEN_ENGLISH
609 
643 #endif
644 #ifdef DOXYGEN_RUSSIAN
645 
683 #endif
684 int DSPL_API ellip_cd_cmplx(complex_t* u, int n, double k, complex_t* y)
685 {
686  double lnd[ELLIP_ITER], t;
687  int i, m;
688  complex_t tmp;
689 
690  if(!u || !y)
691  return ERROR_PTR;
692  if(n<1)
693  return ERROR_SIZE;
694  if(k < 0.0 || k>= 1.0)
695  return ERROR_ELLIP_MODULE;
696 
697  ellip_landen(k,ELLIP_ITER, lnd);
698 
699  for(m = 0; m < n; m++)
700  {
701  RE(tmp) = RE(u[m]) * M_PI * 0.5;
702  IM(tmp) = IM(u[m]) * M_PI * 0.5;
703 
704  cos_cmplx(&tmp, 1, y+m);
705 
706  for(i = ELLIP_ITER-1; i>0; i--)
707  {
708  t = 1.0 / ABSSQR(y[m]);
709 
710  RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i];
711  IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i];
712 
713  t = (1.0 + lnd[i]) / ABSSQR(tmp);
714 
715  RE(y[m]) = RE(tmp) * t;
716  IM(y[m]) = -IM(tmp) * t;
717  }
718  }
719  return RES_OK;
720 }
721 
722 
723 
724 
725 
726 
727 #ifdef DOXYGEN_ENGLISH
728 
797 #endif
798 #ifdef DOXYGEN_RUSSIAN
799 
871 #endif
872 int DSPL_API ellip_landen(double k, int n, double* y)
873 {
874  int i;
875  y[0] = k;
876 
877  if(!y)
878  return ERROR_PTR;
879  if(n < 1)
880  return ERROR_SIZE;
881  if(k < 0.0 || k>= 1.0)
882  return ERROR_ELLIP_MODULE;
883 
884  for(i = 1; i < n; i++)
885  {
886  y[i] = y[i-1] / (1.0 + sqrt(1.0 - y[i-1] * y[i-1]));
887  y[i] *= y[i];
888  }
889 
890  return RES_OK;
891 }
892 
893 
894 
895 
896 
897 #ifdef DOXYGEN_ENGLISH
898 
899 #endif
900 #ifdef DOXYGEN_RUSSIAN
901 
902 #endif
903 int DSPL_API ellip_modulareq(double rp, double rs, int ord, double *k)
904 {
905  double ep, es, ke, kp, t, sn = 0.0;
906  int i, L, r;
907 
908  if(rp < 0 || rp == 0)
909  return ERROR_FILTER_RP;
910  if(rs < 0 || rs == 0)
911  return ERROR_FILTER_RS;
912  if(ord < 1)
913  return ERROR_FILTER_ORD;
914  if(!k)
915  return ERROR_PTR;
916 
917 
918  ep = sqrt(pow(10.0, rp*0.1)-1.0);
919  es = sqrt(pow(10.0, rs*0.1)-1.0);
920 
921  ke = ep/es;
922 
923  ke = sqrt(1.0 - ke*ke);
924 
925  r = ord % 2;
926  L = (ord-r)/2;
927 
928  kp = 1.0;
929  for(i = 0; i < L; i++)
930  {
931  t = (double)(2*i+1) / (double)ord;
932  ellip_sn(&t, 1, ke, &sn);
933  sn*=sn;
934  kp *= sn*sn;
935  }
936 
937  kp *= pow(ke, (double)ord);
938  *k = sqrt(1.0 - kp*kp);
939 
940  return RES_OK;
941 
942 }
943 
944 
945 
946 #ifdef DOXYGEN_ENGLISH
947 
948 #endif
949 #ifdef DOXYGEN_RUSSIAN
950 
951 #endif
952 int DSPL_API ellip_rat(double* w, int n, int ord, double k, double* u)
953 {
954  double t, xi, w2, xi2, k2;
955  int i, m, r, L;
956 
957  if(!u || !w)
958  return ERROR_PTR;
959  if(n<1)
960  return ERROR_SIZE;
961  if(k < 0.0 || k>= 1.0)
962  return ERROR_ELLIP_MODULE;
963 
964  r = ord%2;
965  L = (ord-r)/2;
966 
967  if(r)
968  memcpy(u, w, n*sizeof(double));
969  else
970  {
971  for(m = 0; m < n; m++)
972  {
973  u[m] = 1.0;
974  }
975  }
976 
977  k2 = k*k;
978  for(i = 0; i < L; i++)
979  {
980  t = (double)(2*i+1) / (double)ord;
981  ellip_cd(&t, 1, k, &xi);
982  xi2 = xi*xi;
983  for(m = 0; m < n; m++)
984  {
985  w2 = w[m]*w[m];
986  u[m] *= (w2 - xi2) / (1.0 - w2 * k2 * xi2);
987  u[m] *= (1.0 - k2*xi2) / (1.0 - xi2);
988  }
989  }
990  return RES_OK;
991 }
992 
993 
994 
995 
996 
997 
998 #ifdef DOXYGEN_ENGLISH
999 
1033 #endif
1034 #ifdef DOXYGEN_RUSSIAN
1035 
1073 #endif
1074 int DSPL_API ellip_sn(double* u, int n, double k, double* y)
1075 {
1076  double lnd[ELLIP_ITER];
1077  int i, m;
1078 
1079  if(!u || !y)
1080  return ERROR_PTR;
1081  if(n<1)
1082  return ERROR_SIZE;
1083  if(k < 0.0 || k>= 1.0)
1084  return ERROR_ELLIP_MODULE;
1085 
1086  ellip_landen(k,ELLIP_ITER, lnd);
1087 
1088 
1089  for(m = 0; m < n; m++)
1090  {
1091  y[m] = sin(u[m] * M_PI * 0.5);
1092  for(i = ELLIP_ITER-1; i>0; i--)
1093  {
1094  y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]);
1095  }
1096  }
1097  return RES_OK;
1098 }
1099 
1100 
1101 
1102 
1103 #ifdef DOXYGEN_ENGLISH
1104 
1138 #endif
1139 #ifdef DOXYGEN_RUSSIAN
1140 
1178 #endif
1179 int DSPL_API ellip_sn_cmplx(complex_t* u, int n, double k, complex_t* y)
1180 {
1181  double lnd[ELLIP_ITER], t;
1182  int i, m;
1183  complex_t tmp;
1184 
1185  if(!u || !y)
1186  return ERROR_PTR;
1187  if(n<1)
1188  return ERROR_SIZE;
1189  if(k < 0.0 || k>= 1.0)
1190  return ERROR_ELLIP_MODULE;
1191 
1192  ellip_landen(k,ELLIP_ITER, lnd);
1193 
1194 
1195  for(m = 0; m < n; m++)
1196  {
1197  RE(tmp) = RE(u[m]) * M_PI * 0.5;
1198  IM(tmp) = IM(u[m]) * M_PI * 0.5;
1199 
1200  sin_cmplx(&tmp, 1, y+m);
1201 
1202  for(i = ELLIP_ITER-1; i>0; i--)
1203  {
1204  t = 1.0 / ABSSQR(y[m]);
1205 
1206  RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i];
1207  IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i];
1208 
1209  t = (1.0 + lnd[i]) / ABSSQR(tmp);
1210 
1211  RE(y[m]) = RE(tmp) * t;
1212  IM(y[m]) = -IM(tmp) * t;
1213 
1214  }
1215  }
1216  return RES_OK;
1217 }
1218 
int DSPL_API ellip_cd(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:577
#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:514
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
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:545
#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:553
#define ERROR_FILTER_RP
The filter non-uniformity parameter in the passband is set incorrectly. This parameter is specified i...
Definition: dspl.h:516
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:517
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:1179
int DSPL_API ellip_sn(double *u, int n, double k, double *y)
Jacobi elliptic function of real vector argument.
Definition: ellipj.c:1074
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:497
#define ABSSQR(x)
The macro returns the square of the modulus of a complex number x.
Definition: dspl.h:473
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:684
#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:508
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:417
int DSPL_API ellip_landen(double k, int n, double *y)
Function calculates complete elliptical integral coefficients .
Definition: ellipj.c:872