ARTS: legendre.cc Source File

00001 /* Copyright (C) 2003-2008 Oliver Lemke <olemke@core-dump.info>
00002 
00003    This program is free software; you can redistribute it and/or modify it
00004    under the terms of the GNU General Public License as published by the
00005    Free Software Foundation; either version 2, or (at your option) any
00006    later version.
00007 
00008    This program is distributed in the hope that it will be useful,
00009    but WITHOUT ANY WARRANTY; without even the implied warranty of
00010    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011    GNU General Public License for more details.
00012 
00013    You should have received a copy of the GNU General Public License
00014    along with this program; if not, write to the Free Software
00015    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
00016    USA.
00017 */
00018 
00019 
00020 
00022 //   File description
00024 
00033 #include "legendre.h"
00034 
00035 #ifdef HAVE_SSTREAM
00036 #include <sstream>
00037 #else
00038 #include "sstream.h"
00039 #endif
00040 #include "exceptions.h"
00041 #include <cmath>
00042 
00043 #include "math_funcs.h"
00044 
00045 
00047 
00065 Numeric
00066 legendre_poly (Index l, Index m, Numeric x)
00067 {
00068   Numeric pmm;
00069   Numeric result = 0.;
00070 
00071   if (m < 0 || m > l || abs (x) > 1.0)
00072     {
00073       ostringstream os;
00074       os << "legendre_poly: Condition 0 <= m <= l && -1 < x < 1 failed"
00075         << endl << "  l = " << l << "  m = " << m << "  x = " << x << endl;
00076       throw runtime_error (os.str ());
00077     }
00078 
00079   pmm = 1.0;
00080   if (m > 0)
00081     {
00082       Numeric fact, somx2;
00083 
00084       somx2 = sqrt ((1.0 - x) * (1.0 + x));
00085       fact = 1.0;
00086       for (Index i = 1; i <= m; i++)
00087         {
00088           pmm *= -fact * somx2;
00089           fact += 2.0;
00090         }
00091     }
00092 
00093   if (l == m)
00094     result = pmm;
00095   else
00096     {
00097       Numeric pmmp1;
00098 
00099       pmmp1 = x * (Numeric)(2*m + 1) * pmm;
00100       if (l == (m+1))
00101         result = pmmp1;
00102       else
00103         {
00104           for (Index ll = m+2; ll <= l; ll++)
00105             {
00106               result = (x * (Numeric)(2*ll - 1) * pmmp1 - (Numeric)(ll + m - 1) * pmm) / (Numeric)(ll - m);
00107               pmm = pmmp1;
00108               pmmp1 = result;
00109             }
00110         }
00111     }
00112 
00113   return (result);
00114 }
00115 
00116 
00118 
00136 Numeric
00137 legendre_poly_norm_schmidt (Index l, Index m, Numeric x)
00138 {
00139   Numeric result;
00140 
00141   if (m != 0)
00142 
00143     {
00144       result = ((sqrt (2.0 * fac (l - m) / fac (l + m))
00145                  * legendre_poly (l, m, x)));
00146     }
00147   else
00148     {
00149       result = (legendre_poly (l, m, x));
00150     }
00151 
00152   return(result);
00153 
00154 }
00155 
00156 
00158 
00172 Numeric
00173 legendre_poly_deriv (Index l, Index m, Numeric x)
00174 {
00175   assert (x != 1.);
00176   if (x == 1.)
00177     {
00178       ostringstream os;
00179       os << "legendre_poly_deriv: Condition x != 1 failed"
00180         << endl << "  x = " << x << endl;
00181       throw runtime_error (os.str ());
00182     }
00183   Numeric result;
00184 
00185   if (l == 1)
00186     {
00187       if (m == 0)
00188         {
00189           result = 1;
00190         }
00191       else if (m == 1)
00192         {
00193           result = x / sqrt(1 - x*x);
00194         }
00195       else
00196         {
00197           ostringstream os;
00198           os << "legendre_poly_deriv: "
00199             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00200             << "l = " << l << "  m = " << m << endl;
00201           throw runtime_error (os.str ());
00202         }
00203     }
00204   else if ( m < l )
00205     {
00206       try
00207         {
00208           result = ((Numeric)(l + m) * legendre_poly (l-1, m, x) - 
00209                     (Numeric)l * x * legendre_poly (l, m, x)) /
00210             (1 - x * x);
00211         }
00212       catch (runtime_error e)
00213         {
00214           ostringstream os;
00215           os << e.what () << "legendre_poly_deriv: "
00216             << "Condition m < l failed" << endl
00217             << "l = " << l << "  m = " << m << endl;
00218           throw runtime_error (os.str ());
00219         }
00220     }
00221   else
00222     {
00223       try
00224         {
00225           result = (Numeric)m * x * legendre_poly (l, m, x) / (1 - x * x)
00226             + (Numeric)((l + m) * (l - m + 1)) * legendre_poly (l, m - 1, x)
00227             / sqrt (1 - x * x);
00228         }
00229       catch (runtime_error e)
00230         {
00231           ostringstream os;
00232           os << e.what () << "legendre_poly_norm_schmidt_deriv: "
00233             << "Condition m = l failed" << endl
00234             << "l = " << l << "  m = " << m << endl;
00235           throw runtime_error (os.str ());
00236         }
00237     }
00238  return (result);
00239 }
00240 
00242 
00256 Numeric
00257 legendre_poly_norm_schmidt_deriv (Index l, Index m, Numeric x)
00258 {
00259   assert (x != 1.);
00260   if (x == 1.)
00261     {
00262       ostringstream os;
00263       os << "legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00264         << endl << "  x = " << x << endl;
00265       throw runtime_error (os.str ());
00266     }
00267 
00268   Numeric result;
00269 
00270   if (l == 1)
00271     {
00272       if (m == 0)
00273         {
00274           result = sqrt (2.0 * fac (1 - m) / fac (1 + m));
00275         }
00276       else if (m == 1)
00277         {
00278           result = sqrt (2.0 * fac (1 - m) / fac (1 + m)) * x / sqrt(1 - x*x);
00279         }
00280       else
00281         {
00282           ostringstream os;
00283           os << "legendre_poly_norm_schmidt_deriv: "
00284             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00285             << "l = " << l << "  m = " << m << endl;
00286           throw runtime_error (os.str ());
00287         }
00288 
00289     }
00290   else if ( m < l )
00291     {
00292       try
00293         {
00294           result = ((Numeric)(l + m) * legendre_poly_norm_schmidt (l-1, m, x) - 
00295                     (Numeric)l * x * legendre_poly_norm_schmidt (l, m, x)) /
00296             (1 - x * x);
00297         }
00298       catch (runtime_error e)
00299         {
00300           ostringstream os;
00301           os << e.what () << "legendre_poly_norm_schmidt_deriv: "
00302             << "Condition m < l failed" << endl
00303             << "l = " << l << "  m = " << m << endl;
00304           throw runtime_error (os.str ());
00305         }
00306     }
00307   else
00308     {
00309       try
00310         {
00311           result = (Numeric)m * x * legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
00312             + (Numeric)((l + m) * (l - m + 1)) * legendre_poly_norm_schmidt (l, m - 1, x)
00313             / sqrt (1 - x * x);
00314         }
00315       catch (runtime_error e)
00316         {
00317           ostringstream os;
00318           os << e.what () << "legendre_poly_norm_schmidt_deriv: "
00319             << "Condition m = l failed" << endl
00320             << "l = " << l << "  m = " << m << endl;
00321           throw runtime_error (os.str ());
00322         }
00323     }
00324 
00325   return (result);
00326 }
00327 
00329 
00348 Numeric
00349 g_legendre_poly (Index l, Index m, Numeric x)
00350 {
00351   Numeric pmm;
00352   Numeric result = 0.;
00353 
00354   if (m < 0 || m > l || abs (x) > 1.0)
00355     {
00356       ostringstream os;
00357       os << "g_legendre_poly: Condition 0 <= m <= l && -1 < x < 1 failed"
00358         << endl << "  l = " << l << "  m = " << m << "  x = " << x << endl;
00359       throw runtime_error (os.str ());
00360     }
00361 
00362   pmm = 1.0;
00363   if (m > 0)
00364     {
00365       Numeric fact, somx2;
00366 
00367       somx2 = sqrt ((1.0 - x) * (1.0 + x));
00368       fact = 1.0;
00369       for (Index i = 1; i <= m; i++)
00370         {
00371           pmm *= fact * somx2;
00372           fact += 2.0;
00373         }
00374     }
00375 
00376   if (l == m)
00377     result = pmm;
00378   else
00379     {
00380       Numeric pmmp1;
00381 
00382       pmmp1 = x * (Numeric)(2*m + 1) * pmm;
00383       if (l == (m+1))
00384         result = pmmp1;
00385       else
00386         {
00387           for (Index ll = m+2; ll <= l; ll++)
00388             {
00389               result = (x * (Numeric)(2*ll - 1) * pmmp1 - (Numeric)(ll + m - 1) * pmm) / (Numeric)(ll - m);
00390               pmm = pmmp1;
00391               pmmp1 = result;
00392             }
00393         }
00394     }
00395 
00396   return (result);
00397 }
00398 
00399 
00401 
00420 Numeric
00421 g_legendre_poly_norm_schmidt (Index l, Index m, Numeric x)
00422 {
00423   Numeric result;
00424 
00425   if (m != 0)
00426 
00427     {
00428       result = ((sqrt (2.0 * fac (l - m) / fac (l + m))
00429                  * g_legendre_poly (l, m, x)));
00430     }
00431   else
00432     {
00433       result = (g_legendre_poly (l, m, x));
00434     }
00435 
00436   return(result);
00437 
00438 }
00439 
00440 
00442 
00457 Numeric
00458 g_legendre_poly_deriv (Index l, Index m, Numeric x)
00459 {
00460   assert (x != 1.);
00461   if (x == 1.)
00462     {
00463       ostringstream os;
00464       os << "g_legendre_poly_deriv: Condition x != 1 failed"
00465         << endl << "  x = " << x << endl;
00466       throw runtime_error (os.str ());
00467     }
00468   Numeric result;
00469 
00470   if (l == 1)
00471     {
00472       if (m == 0)
00473         {
00474           result = 1;
00475         }
00476       else if (m == 1)
00477         {
00478           result = x / sqrt(1 - x*x);
00479         }
00480       else
00481         {
00482           ostringstream os;
00483           os << "g_legendre_poly_deriv: "
00484             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00485             << "l = " << l << "  m = " << m << endl;
00486           throw runtime_error (os.str ());
00487         }
00488     }
00489   else if ( m < l )
00490     {
00491       try
00492         {
00493           result = ((Numeric)(l + m) * g_legendre_poly (l-1, m, x) - 
00494                     (Numeric)l * x * g_legendre_poly (l, m, x)) /
00495             (1 - x * x);
00496         }
00497       catch (runtime_error e)
00498         {
00499           ostringstream os;
00500           os << e.what () << "g_legendre_poly_deriv: "
00501             << "Condition m < l failed" << endl
00502             << "l = " << l << "  m = " << m << endl;
00503           throw runtime_error (os.str ());
00504         }
00505     }
00506   else
00507     {
00508       try
00509         {
00510           result = - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x) + 
00511             (Numeric)((l + m) * (l - m + 1)) * g_legendre_poly (l, m - 1, x) / 
00512             sqrt (1 - x * x);
00513         }
00514       catch (runtime_error e)
00515         {
00516           ostringstream os;
00517           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00518             << "Condition m = l failed" << endl
00519             << "l = " << l << "  m = " << m << endl;
00520           throw runtime_error (os.str ());
00521         }
00522     }
00523  return (result);
00524 }
00525 
00527 
00542 Numeric
00543 g_legendre_poly_norm_schmidt_deriv (Index l, Index m, Numeric x)
00544 {
00545   assert (x != 1.);
00546   if (x == 1.)
00547     {
00548       ostringstream os;
00549       os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00550         << endl << "  x = " << x << endl;
00551       throw runtime_error (os.str ());
00552     }
00553 
00554   Numeric result;
00555 
00556   if (l == 1)
00557     {
00558       if (m == 0)
00559         {
00560           result = 1;
00561         }
00562       else if (m == 1)
00563         {
00564           result = x / sqrt(1 - x * x);
00565         }
00566       else
00567         {
00568           ostringstream os;
00569           os << "g_legendre_poly_norm_schmidt_deriv: "
00570             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00571             << "l = " << l << "  m = " << m << endl;
00572           throw runtime_error (os.str ());
00573         }
00574 
00575     }
00576   else if ( m < l )
00577     {
00578       try
00579         {
00580           result = sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00581             ((Numeric)(l + m) * g_legendre_poly(l-1, m, x) - 
00582              (Numeric)l * x * g_legendre_poly (l, m, x)) / (1 - x * x);
00583         }
00584       catch (runtime_error e)
00585         {
00586           ostringstream os;
00587           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00588             << "Condition m < l failed" << endl
00589             << "l = " << l << "  m = " << m << endl;
00590           throw runtime_error (os.str ());
00591         }
00592     }
00593   else
00594     {
00595       try
00596         {
00597           result = sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00598             ( - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x) + 
00599               (Numeric)((l + m) * (l - m + 1)) * g_legendre_poly (l, m - 1, x) / 
00600               sqrt (1 - x * x));
00601         }
00602       catch (runtime_error e)
00603         {
00604           ostringstream os;
00605           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00606             << "Condition m = l failed" << endl
00607             << "l = " << l << "  m = " << m << endl;
00608           throw runtime_error (os.str ());
00609         }
00610     }
00611 
00612   return (result);
00613 }
00614 
00616 
00631 Numeric
00632 g_legendre_poly_norm_schmidt_deriv1 (Index l, Index m, Numeric x)
00633 {
00634   assert (x != 1.);
00635   if (x == 1.)
00636     {
00637       ostringstream os;
00638       os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00639         << endl << "  x = " << x << endl;
00640       throw runtime_error (os.str ());
00641     }
00642 
00643   Numeric result;
00644 
00645   if (l == 1)
00646     {
00647       if (m == 0)
00648         {
00649           result = 1;
00650         }
00651       else if (m == 1)
00652         {
00653           result = x / sqrt(1 - x * x);
00654         }
00655       else
00656         {
00657           ostringstream os;
00658           os << "g_legendre_poly_norm_schmidt_deriv: "
00659             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00660             << "l = " << l << "  m = " << m << endl;
00661           throw runtime_error (os.str ());
00662         }
00663 
00664     }
00665   else if ( m <= l - 1 )
00666     {
00667       try
00668         {
00669           // result = - m * x * g_legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
00670           //  +  sqrt((double)(l + m + 1 / l - m - 1)) * g_legendre_poly_norm_schmidt (l, m + 1, x)
00671           //  / sqrt (1 - x * x);
00672 
00673           result = sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00674             ( - (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x)  + 
00675               g_legendre_poly (l, m + 1, x)
00676               / sqrt (1 - x * x));
00677         }
00678       catch (runtime_error e)
00679         {
00680           ostringstream os;
00681           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00682             << "Condition m <= l - 1 failed" << endl
00683             << "l = " << l << "  m = " << m << endl;
00684           throw runtime_error (os.str ());
00685         }
00686     }
00687   else
00688     {
00689       try
00690         {
00691           result = - sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00692             (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
00693         }
00694       catch (runtime_error e)
00695         {
00696           ostringstream os;
00697           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00698             << "Condition m = l failed" << endl
00699             << "l = " << l << "  m = " << m << endl;
00700           throw runtime_error (os.str ());
00701         }
00702     }
00703 
00704   return (result);
00705 }
00706 
00708 
00723 Numeric
00724 g_legendre_poly_norm_schmidt_deriv2 (Index l, Index m, Numeric x)
00725 {
00726   assert (x != 1.);
00727   if (x == 1.)
00728     {
00729       ostringstream os;
00730       os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00731         << endl << "  x = " << x << endl;
00732       throw runtime_error (os.str ());
00733     }
00734 
00735   Numeric result;
00736 
00737   if (l == 1)
00738     {
00739       if (m == 0)
00740         {
00741           result = 1;
00742         }
00743       else if (m == 1)
00744         {
00745           result = x / sqrt(1 - x * x);
00746         }
00747       else
00748         {
00749           ostringstream os;
00750           os << "g_legendre_poly_norm_schmidt_deriv: "
00751             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00752             << "l = " << l << "  m = " << m << endl;
00753           throw runtime_error (os.str ());
00754         }
00755 
00756     }
00757   else if ( m < l )
00758     {
00759       try
00760         {
00761           result = - sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00762             ((Numeric)(l + m) * g_legendre_poly (l-1, m, x) - 
00763              (Numeric)l * x * g_legendre_poly (l, m, x)) /
00764             (1 - x * x);
00765         }
00766       catch (runtime_error e)
00767         {
00768           ostringstream os;
00769           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00770             << "Condition m < l failed" << endl
00771             << "l = " << l << "  m = " << m << endl;
00772           throw runtime_error (os.str ());
00773         }
00774     }
00775   else
00776     {
00777       try
00778         {
00779           result = - (Numeric)m * x * g_legendre_poly_norm_schmidt (l, m, x) / (1 - x * x)
00780             ;
00781         }
00782       catch (runtime_error e)
00783         {
00784           ostringstream os;
00785           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00786             << "Condition m = l failed" << endl
00787             << "l = " << l << "  m = " << m << endl;
00788           throw runtime_error (os.str ());
00789         }
00790     }
00791 
00792   return (result);
00793 }
00794 
00796 
00811 Numeric
00812 g_legendre_poly_norm_schmidt_deriv3 (Index l, Index m, Numeric x)
00813 {
00814   assert (x != 1.);
00815   if (x == 1.)
00816     {
00817       ostringstream os;
00818       os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00819         << endl << "  x = " << x << endl;
00820       throw runtime_error (os.str ());
00821     }
00822 
00823   Numeric result;
00824 
00825   if (l == 1)
00826     {
00827       if (m == 0)
00828         {
00829           result = 1;
00830         }
00831       else if (m == 1)
00832         {
00833           result = x / sqrt(1 - x * x);
00834         }
00835       else
00836         {
00837           ostringstream os;
00838           os << "g_legendre_poly_norm_schmidt_deriv: "
00839             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00840             << "l = " << l << "  m = " << m << endl;
00841           throw runtime_error (os.str ());
00842         }
00843 
00844     }
00845   else if ( m < l )
00846     {
00847       try
00848         {
00849           result = sqrt(2.0 * fac (l - m) / fac (l + m)) * 
00850             ((Numeric)l * g_legendre_poly (l - 1, m, x) + 
00851              (Numeric)(m - l) * x * g_legendre_poly (l, m, x)) /
00852             (1 - x * x);
00853         }
00854       catch (runtime_error e)
00855         {
00856           ostringstream os;
00857           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00858             << "Condition m < l failed" << endl
00859             << "l = " << l << "  m = " << m << endl;
00860           throw runtime_error (os.str ());
00861         }
00862     }
00863   else
00864     {
00865       try
00866         {
00867           result = - sqrt(2.0 * fac (l - m) / fac (l + m)) * 
00868             (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
00869         }
00870       catch (runtime_error e)
00871         {
00872           ostringstream os;
00873           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00874             << "Condition m = l failed" << endl
00875             << "l = " << l << "  m = " << m << endl;
00876           throw runtime_error (os.str ());
00877         }
00878     }
00879 
00880   return (result);
00881 }
00882 
00884 
00899 Numeric
00900 g_legendre_poly_norm_schmidt_deriv4 (Index l, Index m, Numeric x)
00901 {
00902   assert (x != 1.);
00903   if (x == 1.)
00904     {
00905       ostringstream os;
00906       os << "g_legendre_poly_norm_schmidt_deriv: Condition x != 1 failed"
00907         << endl << "  x = " << x << endl;
00908       throw runtime_error (os.str ());
00909     }
00910 
00911   Numeric result;
00912 
00913   if (l == 1)
00914     {
00915       if (m == 0)
00916         {
00917           result = 1;
00918         }
00919       else if (m == 1)
00920         {
00921           result = x / sqrt(1 - x * x);
00922         }
00923       else
00924         {
00925           ostringstream os;
00926           os << "g_legendre_poly_norm_schmidt_deriv: "
00927             << "Condition l == 1 && (m == 0 || m == 1) failed" << endl
00928             << "l = " << l << "  m = " << m << endl;
00929           throw runtime_error (os.str ());
00930         }
00931 
00932     }
00933   else if ( m < l )
00934     {
00935       try
00936         {
00937           result = sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00938             ((Numeric)((l + m) * (l + 1)) * g_legendre_poly (l - 1, m, x)  - 
00939              (Numeric)((l + 2 * m) * (l - m + 1)) * g_legendre_poly (l + 1, m, x)
00940              / ((Numeric)(2 * l + 1 ) * (1 - x * x)));
00941         }
00942       catch (runtime_error e)
00943         {
00944           ostringstream os;
00945           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00946             << "Condition m < l failed" << endl
00947             << "l = " << l << "  m = " << m << endl;
00948           throw runtime_error (os.str ());
00949         }
00950     }
00951   else
00952     {
00953       try
00954         {
00955           result = - sqrt (2.0 * fac (l - m) / fac (l + m)) * 
00956             (Numeric)m * x * g_legendre_poly (l, m, x) / (1 - x * x);
00957         }
00958       catch (runtime_error e)
00959         {
00960           ostringstream os;
00961           os << e.what () << "g_legendre_poly_norm_schmidt_deriv: "
00962             << "Condition m = l failed" << endl
00963             << "l = " << l << "  m = " << m << endl;
00964           throw runtime_error (os.str ());
00965         }
00966     }
00967 
00968   return (result);
00969 }