ARTS  2.0.49
legendre.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2008 Oliver Lemke <olemke@core-dump.info>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA.
17 */
18 
19 
20 
22 // File description
24 
33 #include "legendre.h"
34 
35 #ifdef HAVE_SSTREAM
36 #include <sstream>
37 #else
38 #include "sstream.h"
39 #endif
40 #include "exceptions.h"
41 #include <cmath>
42 
43 #include "math_funcs.h"
44 
45 
47 
65 Numeric
67 {
68  Numeric pmm;
69  Numeric result = 0.;
70 
71  if (m < 0 || m > l || abs (x) > 1.0)
72  {
73  ostringstream os;
74  os << "legendre_poly: Condition 0 <= m <= l && -1 < x < 1 failed"
75  << endl << " l = " << l << " m = " << m << " x = " << x << endl;
76  throw runtime_error (os.str ());
77  }
78 
79  pmm = 1.0;
80  if (m > 0)
81  {
82  Numeric fact, somx2;
83 
84  somx2 = sqrt ((1.0 - x) * (1.0 + x));
85  fact = 1.0;
86  for (Index i = 1; i <= m; i++)
87  {
88  pmm *= -fact * somx2;
89  fact += 2.0;
90  }
91  }
92 
93  if (l == m)
94  result = pmm;
95  else
96  {
97  Numeric pmmp1;
98 
99  pmmp1 = x * (Numeric)(2*m + 1) * pmm;
100  if (l == (m+1))
101  result = pmmp1;
102  else
103  {
104  for (Index ll = m+2; ll <= l; ll++)
105  {
106  result = (x * (Numeric)(2*ll - 1) * pmmp1 - (Numeric)(ll + m - 1) * pmm) / (Numeric)(ll - m);
107  pmm = pmmp1;
108  pmmp1 = result;
109  }
110  }
111  }
112 
113  return (result);
114 }
115 
116 
118 
136 Numeric
138 {
139  Numeric result;
140 
141  if (m != 0)
142 
143  {
144  result = ((sqrt (2.0 * fac (l - m) / fac (l + m))
145  * legendre_poly (l, m, x)));
146  }
147  else
148  {
149  result = (legendre_poly (l, m, x));
150  }
151 
152  return(result);
153 
154 }
155 
156 
158 
172 Numeric
174 {
175  assert (x != 1.);
176  if (x == 1.)
177  {
178  ostringstream os;
179  os << "legendre_poly_deriv: Condition x != 1 failed"
180  << endl << " x = " << x << endl;
181  throw runtime_error (os.str ());
182  }
183  Numeric result;
184 
185  if (l == 1)
186  {
187  if (m == 0)
188  {
189  result = 1;
190  }
191  else if (m == 1)
192  {
193  result = x / sqrt(1 - x*x);
194  }
195  else
196  {
197  ostringstream os;
198  os << "legendre_poly_deriv: "
199  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
200  << "l = " << l << " m = " << m << endl;
201  throw runtime_error (os.str ());
202  }
203  }
204  else if ( m < l )
205  {
206  try
207  {
208  result = ((Numeric)(l + m) * legendre_poly (l-1, m, x) -
209  (Numeric)l * x * legendre_poly (l, m, x)) /
210  (1 - x * x);
211  }
212  catch (runtime_error e)
213  {
214  ostringstream os;
215  os << e.what () << "legendre_poly_deriv: "
216  << "Condition m < l failed" << endl
217  << "l = " << l << " m = " << m << endl;
218  throw runtime_error (os.str ());
219  }
220  }
221  else
222  {
223  try
224  {
225  result = (Numeric)m * x * legendre_poly (l, m, x) / (1 - x * x)
226  + (Numeric)((l + m) * (l - m + 1)) * legendre_poly (l, m - 1, x)
227  / sqrt (1 - x * x);
228  }
229  catch (runtime_error e)
230  {
231  ostringstream os;
232  os << e.what () << "legendre_poly_norm_schmidt_deriv: "
233  << "Condition m = l failed" << endl
234  << "l = " << l << " m = " << m << endl;
235  throw runtime_error (os.str ());
236  }
237  }
238  return (result);
239 }
240 
242 
256 Numeric
258 {
259  assert (x != 1.);
260  if (x == 1.)
261  {
262  ostringstream os;
263  os << "legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
264  << endl << " x = " << x << endl;
265  throw runtime_error (os.str ());
266  }
267 
268  Numeric result;
269 
270  if (l == 1)
271  {
272  if (m == 0)
273  {
274  result = sqrt (2.0 * fac (1 - m) / fac (1 + m));
275  }
276  else if (m == 1)
277  {
278  result = sqrt (2.0 * fac (1 - m) / fac (1 + m)) * x / sqrt(1 - x*x);
279  }
280  else
281  {
282  ostringstream os;
283  os << "legendre_poly_norm_schmidt_deriv: "
284  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
285  << "l = " << l << " m = " << m << endl;
286  throw runtime_error (os.str ());
287  }
288 
289  }
290  else if ( m < l )
291  {
292  try
293  {
294  result = ((Numeric)(l + m) * legendre_poly_norm_schmidt (l-1, m, x) -
295  (Numeric)l * x * legendre_poly_norm_schmidt (l, m, x)) /
296  (1 - x * x);
297  }
298  catch (runtime_error e)
299  {
300  ostringstream os;
301  os << e.what () << "legendre_poly_norm_schmidt_deriv: "
302  << "Condition m < l failed" << endl
303  << "l = " << l << " m = " << m << endl;
304  throw runtime_error (os.str ());
305  }
306  }
307  else
308  {
309  try
310  {
311  result = (Numeric)m * x * legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
312  + (Numeric)((l + m) * (l - m + 1)) * legendre_poly_norm_schmidt (l, m - 1, x)
313  / sqrt (1 - x * x);
314  }
315  catch (runtime_error e)
316  {
317  ostringstream os;
318  os << e.what () << "legendre_poly_norm_schmidt_deriv: "
319  << "Condition m = l failed" << endl
320  << "l = " << l << " m = " << m << endl;
321  throw runtime_error (os.str ());
322  }
323  }
324 
325  return (result);
326 }
327 
329 
348 Numeric
350 {
351  Numeric pmm;
352  Numeric result = 0.;
353 
354  if (m < 0 || m > l || abs (x) > 1.0)
355  {
356  ostringstream os;
357  os << "g_legendre_poly: Condition 0 <= m <= l && -1 < x < 1 failed"
358  << endl << " l = " << l << " m = " << m << " x = " << x << endl;
359  throw runtime_error (os.str ());
360  }
361 
362  pmm = 1.0;
363  if (m > 0)
364  {
365  Numeric fact, somx2;
366 
367  somx2 = sqrt ((1.0 - x) * (1.0 + x));
368  fact = 1.0;
369  for (Index i = 1; i <= m; i++)
370  {
371  pmm *= fact * somx2;
372  fact += 2.0;
373  }
374  }
375 
376  if (l == m)
377  result = pmm;
378  else
379  {
380  Numeric pmmp1;
381 
382  pmmp1 = x * (Numeric)(2*m + 1) * pmm;
383  if (l == (m+1))
384  result = pmmp1;
385  else
386  {
387  for (Index ll = m+2; ll <= l; ll++)
388  {
389  result = (x * (Numeric)(2*ll - 1) * pmmp1 - (Numeric)(ll + m - 1) * pmm) / (Numeric)(ll - m);
390  pmm = pmmp1;
391  pmmp1 = result;
392  }
393  }
394  }
395 
396  return (result);
397 }
398 
399 
401 
420 Numeric
422 {
423  Numeric result;
424 
425  if (m != 0)
426 
427  {
428  result = ((sqrt (2.0 * fac (l - m) / fac (l + m))
429  * g_legendre_poly (l, m, x)));
430  }
431  else
432  {
433  result = (g_legendre_poly (l, m, x));
434  }
435 
436  return(result);
437 
438 }
439 
440 
442 
457 Numeric
459 {
460  assert (x != 1.);
461  if (x == 1.)
462  {
463  ostringstream os;
464  os << "g_legendre_poly_deriv: Condition x != 1 failed"
465  << endl << " x = " << x << endl;
466  throw runtime_error (os.str ());
467  }
468  Numeric result;
469 
470  if (l == 1)
471  {
472  if (m == 0)
473  {
474  result = 1;
475  }
476  else if (m == 1)
477  {
478  result = x / sqrt(1 - x*x);
479  }
480  else
481  {
482  ostringstream os;
483  os << "g_legendre_poly_deriv: "
484  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
485  << "l = " << l << " m = " << m << endl;
486  throw runtime_error (os.str ());
487  }
488  }
489  else if ( m < l )
490  {
491  try
492  {
493  result = ((Numeric)(l + m) * g_legendre_poly (l-1, m, x) -
494  (Numeric)l * x * g_legendre_poly (l, m, x)) /
495  (1 - x * x);
496  }
497  catch (runtime_error e)
498  {
499  ostringstream os;
500  os << e.what () << "g_legendre_poly_deriv: "
501  << "Condition m < l failed" << endl
502  << "l = " << l << " m = " << m << endl;
503  throw runtime_error (os.str ());
504  }
505  }
506  else
507  {
508  try
509  {
510  result = - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x) +
511  (Numeric)((l + m) * (l - m + 1)) * g_legendre_poly (l, m - 1, x) /
512  sqrt (1 - x * x);
513  }
514  catch (runtime_error e)
515  {
516  ostringstream os;
517  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
518  << "Condition m = l failed" << endl
519  << "l = " << l << " m = " << m << endl;
520  throw runtime_error (os.str ());
521  }
522  }
523  return (result);
524 }
525 
527 
542 Numeric
544 {
545  assert (x != 1.);
546  if (x == 1.)
547  {
548  ostringstream os;
549  os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
550  << endl << " x = " << x << endl;
551  throw runtime_error (os.str ());
552  }
553 
554  Numeric result;
555 
556  if (l == 1)
557  {
558  if (m == 0)
559  {
560  result = 1;
561  }
562  else if (m == 1)
563  {
564  result = x / sqrt(1 - x * x);
565  }
566  else
567  {
568  ostringstream os;
569  os << "g_legendre_poly_norm_schmidt_deriv: "
570  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
571  << "l = " << l << " m = " << m << endl;
572  throw runtime_error (os.str ());
573  }
574 
575  }
576  else if ( m < l )
577  {
578  try
579  {
580  result = sqrt (2.0 * fac (l - m) / fac (l + m)) *
581  ((Numeric)(l + m) * g_legendre_poly(l-1, m, x) -
582  (Numeric)l * x * g_legendre_poly (l, m, x)) / (1 - x * x);
583  }
584  catch (runtime_error e)
585  {
586  ostringstream os;
587  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
588  << "Condition m < l failed" << endl
589  << "l = " << l << " m = " << m << endl;
590  throw runtime_error (os.str ());
591  }
592  }
593  else
594  {
595  try
596  {
597  result = sqrt (2.0 * fac (l - m) / fac (l + m)) *
598  ( - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x) +
599  (Numeric)((l + m) * (l - m + 1)) * g_legendre_poly (l, m - 1, x) /
600  sqrt (1 - x * x));
601  }
602  catch (runtime_error e)
603  {
604  ostringstream os;
605  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
606  << "Condition m = l failed" << endl
607  << "l = " << l << " m = " << m << endl;
608  throw runtime_error (os.str ());
609  }
610  }
611 
612  return (result);
613 }
614 
616 
631 Numeric
633 {
634  assert (x != 1.);
635  if (x == 1.)
636  {
637  ostringstream os;
638  os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
639  << endl << " x = " << x << endl;
640  throw runtime_error (os.str ());
641  }
642 
643  Numeric result;
644 
645  if (l == 1)
646  {
647  if (m == 0)
648  {
649  result = 1;
650  }
651  else if (m == 1)
652  {
653  result = x / sqrt(1 - x * x);
654  }
655  else
656  {
657  ostringstream os;
658  os << "g_legendre_poly_norm_schmidt_deriv: "
659  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
660  << "l = " << l << " m = " << m << endl;
661  throw runtime_error (os.str ());
662  }
663 
664  }
665  else if ( m <= l - 1 )
666  {
667  try
668  {
669  // result = - m * x * g_legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
670  // + sqrt((double)(l + m + 1 / l - m - 1)) * g_legendre_poly_norm_schmidt (l, m + 1, x)
671  // / sqrt (1 - x * x);
672 
673  result = sqrt (2.0 * fac (l - m) / fac (l + m)) *
674  ( - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x) +
675  g_legendre_poly (l, m + 1, x)
676  / sqrt (1 - x * x));
677  }
678  catch (runtime_error e)
679  {
680  ostringstream os;
681  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
682  << "Condition m <= l - 1 failed" << endl
683  << "l = " << l << " m = " << m << endl;
684  throw runtime_error (os.str ());
685  }
686  }
687  else
688  {
689  try
690  {
691  result = - sqrt (2.0 * fac (l - m) / fac (l + m)) *
692  (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
693  }
694  catch (runtime_error e)
695  {
696  ostringstream os;
697  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
698  << "Condition m = l failed" << endl
699  << "l = " << l << " m = " << m << endl;
700  throw runtime_error (os.str ());
701  }
702  }
703 
704  return (result);
705 }
706 
708 
723 Numeric
725 {
726  assert (x != 1.);
727  if (x == 1.)
728  {
729  ostringstream os;
730  os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
731  << endl << " x = " << x << endl;
732  throw runtime_error (os.str ());
733  }
734 
735  Numeric result;
736 
737  if (l == 1)
738  {
739  if (m == 0)
740  {
741  result = 1;
742  }
743  else if (m == 1)
744  {
745  result = x / sqrt(1 - x * x);
746  }
747  else
748  {
749  ostringstream os;
750  os << "g_legendre_poly_norm_schmidt_deriv: "
751  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
752  << "l = " << l << " m = " << m << endl;
753  throw runtime_error (os.str ());
754  }
755 
756  }
757  else if ( m < l )
758  {
759  try
760  {
761  result = - sqrt (2.0 * fac (l - m) / fac (l + m)) *
762  ((Numeric)(l + m) * g_legendre_poly (l-1, m, x) -
763  (Numeric)l * x * g_legendre_poly (l, m, x)) /
764  (1 - x * x);
765  }
766  catch (runtime_error e)
767  {
768  ostringstream os;
769  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
770  << "Condition m < l failed" << endl
771  << "l = " << l << " m = " << m << endl;
772  throw runtime_error (os.str ());
773  }
774  }
775  else
776  {
777  try
778  {
779  result = - (Numeric)m * x * g_legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
780  ;
781  }
782  catch (runtime_error e)
783  {
784  ostringstream os;
785  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
786  << "Condition m = l failed" << endl
787  << "l = " << l << " m = " << m << endl;
788  throw runtime_error (os.str ());
789  }
790  }
791 
792  return (result);
793 }
794 
796 
811 Numeric
813 {
814  assert (x != 1.);
815  if (x == 1.)
816  {
817  ostringstream os;
818  os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
819  << endl << " x = " << x << endl;
820  throw runtime_error (os.str ());
821  }
822 
823  Numeric result;
824 
825  if (l == 1)
826  {
827  if (m == 0)
828  {
829  result = 1;
830  }
831  else if (m == 1)
832  {
833  result = x / sqrt(1 - x * x);
834  }
835  else
836  {
837  ostringstream os;
838  os << "g_legendre_poly_norm_schmidt_deriv: "
839  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
840  << "l = " << l << " m = " << m << endl;
841  throw runtime_error (os.str ());
842  }
843 
844  }
845  else if ( m < l )
846  {
847  try
848  {
849  result = sqrt(2.0 * fac (l - m) / fac (l + m)) *
850  ((Numeric)l * g_legendre_poly (l - 1, m, x) +
851  (Numeric)(m - l) * x * g_legendre_poly (l, m, x)) /
852  (1 - x * x);
853  }
854  catch (runtime_error e)
855  {
856  ostringstream os;
857  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
858  << "Condition m < l failed" << endl
859  << "l = " << l << " m = " << m << endl;
860  throw runtime_error (os.str ());
861  }
862  }
863  else
864  {
865  try
866  {
867  result = - sqrt(2.0 * fac (l - m) / fac (l + m)) *
868  (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
869  }
870  catch (runtime_error e)
871  {
872  ostringstream os;
873  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
874  << "Condition m = l failed" << endl
875  << "l = " << l << " m = " << m << endl;
876  throw runtime_error (os.str ());
877  }
878  }
879 
880  return (result);
881 }
882 
884 
899 Numeric
901 {
902  assert (x != 1.);
903  if (x == 1.)
904  {
905  ostringstream os;
906  os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
907  << endl << " x = " << x << endl;
908  throw runtime_error (os.str ());
909  }
910 
911  Numeric result;
912 
913  if (l == 1)
914  {
915  if (m == 0)
916  {
917  result = 1;
918  }
919  else if (m == 1)
920  {
921  result = x / sqrt(1 - x * x);
922  }
923  else
924  {
925  ostringstream os;
926  os << "g_legendre_poly_norm_schmidt_deriv: "
927  << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
928  << "l = " << l << " m = " << m << endl;
929  throw runtime_error (os.str ());
930  }
931 
932  }
933  else if ( m < l )
934  {
935  try
936  {
937  result = sqrt (2.0 * fac (l - m) / fac (l + m)) *
938  ((Numeric)((l + m) * (l + 1)) * g_legendre_poly (l - 1, m, x) -
939  (Numeric)((l + 2 * m) * (l - m + 1)) * g_legendre_poly (l + 1, m, x)
940  / ((Numeric)(2 * l + 1 ) * (1 - x * x)));
941  }
942  catch (runtime_error e)
943  {
944  ostringstream os;
945  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
946  << "Condition m < l failed" << endl
947  << "l = " << l << " m = " << m << endl;
948  throw runtime_error (os.str ());
949  }
950  }
951  else
952  {
953  try
954  {
955  result = - sqrt (2.0 * fac (l - m) / fac (l + m)) *
956  (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
957  }
958  catch (runtime_error e)
959  {
960  ostringstream os;
961  os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
962  << "Condition m = l failed" << endl
963  << "l = " << l << " m = " << m << endl;
964  throw runtime_error (os.str ());
965  }
966  }
967 
968  return (result);
969 }
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:68
exceptions.h
The declarations of all the exception classes.
g_legendre_poly_deriv
Numeric g_legendre_poly_deriv(Index l, Index m, Numeric x)
g_legendre_poly_deriv
Definition: legendre.cc:458
legendre.h
Contains the code to calculate Legendre polynomials.
g_legendre_poly_norm_schmidt_deriv3
Numeric g_legendre_poly_norm_schmidt_deriv3(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv3
Definition: legendre.cc:812
legendre_poly_norm_schmidt
Numeric legendre_poly_norm_schmidt(Index l, Index m, Numeric x)
legendre_poly_norm_schmidt
Definition: legendre.cc:137
g_legendre_poly_norm_schmidt_deriv1
Numeric g_legendre_poly_norm_schmidt_deriv1(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv1
Definition: legendre.cc:632
legendre_poly_norm_schmidt_deriv
Numeric legendre_poly_norm_schmidt_deriv(Index l, Index m, Numeric x)
legendre_poly_norm_schmidt_deriv
Definition: legendre.cc:257
abs
#define abs(x)
Definition: continua.cc:13094
legendre_poly_deriv
Numeric legendre_poly_deriv(Index l, Index m, Numeric x)
legendre_poly_deriv
Definition: legendre.cc:173
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
g_legendre_poly_norm_schmidt
Numeric g_legendre_poly_norm_schmidt(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt
Definition: legendre.cc:421
math_funcs.h
legendre_poly
Numeric legendre_poly(Index l, Index m, Numeric x)
legendre_poly
Definition: legendre.cc:66
g_legendre_poly_norm_schmidt_deriv4
Numeric g_legendre_poly_norm_schmidt_deriv4(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv4
Definition: legendre.cc:900
g_legendre_poly_norm_schmidt_deriv
Numeric g_legendre_poly_norm_schmidt_deriv(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv
Definition: legendre.cc:543
g_legendre_poly_norm_schmidt_deriv2
Numeric g_legendre_poly_norm_schmidt_deriv2(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv2
Definition: legendre.cc:724
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
g_legendre_poly
Numeric g_legendre_poly(Index l, Index m, Numeric x)
g_legendre_poly
Definition: legendre.cc:349
ll
#define ll
Definition: continua.cc:14339