ARTS  2.2.66
wigner_functions.cc
Go to the documentation of this file.
1 /* Copyright (C) 2012
2  * Richard Larsson <ric.larsson@gmail.com>
3  *
4  * This program is free software; you can redistribute it and/or modify it
5  * under the terms of the GNU General Public License as published by the
6  * Free Software Foundation; either version 2, or (at your option) any
7  * later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  * USA. */
18 
19 #include "wigner_functions.h"
20 #include <stdexcept>
21 #include <sstream>
22 
28 Numeric wigner3j(const Rational j1,const Rational j2,const Rational j3,
29  const Rational m1,const Rational m2,const Rational m3)
30 {
31  //std::cout<<std::endl<<"Wigner3j (" << j1 <<" "<<j2 <<" "<<j3<<"; " << m1 <<" " << m2 << " "<<m3<<")"<<std::endl;
32 
33  Rational J = j1 + j2 + j3;
34 
35  {// TEST Area
36  if( !( ( m1 + m2 + m3 ).toIndex() == 0 ) )
37  return 0;
38 
39  if( !( triangular_inequality( j1, j2, j3 ) ) )
40  return 0;
41 
42  if( !( ( J ).Denom() == 1 ) )
43  return 0;
44 
45  if( !( abs( m1 ) <= j1 && abs( m2 ) <= j2 && abs( m3 ) <= j3 ) )
46  return 0;
47 
48  if(m1.toIndex() ==0 && m2.toIndex() == 0 && m3.toIndex() ==0)
49  {
50  if( !( (( J ).toIndex() % 2) == 0 ) )
51  return 0;
52  else
53  {
54  return ((((J/2)%2)==0)?(1.):(-1.)) * sqrt( factorials(
55  MakeArray<Index>( (J-2*j1).toIndex(),
56  (J-2*j2).toIndex(),
57  (J-2*j3).toIndex(),
58  (J/2).toIndex(),
59  (J/2).toIndex()),
60  MakeArray<Index>( (J+1).toIndex(),
61  (J/2-j1).toIndex(),
62  (J/2-j1).toIndex(),
63  (J/2-j2).toIndex(),
64  (J/2-j2).toIndex(),
65  (J/2-j3).toIndex(),
66  (J/2-j3).toIndex())));
67  }
68  }
69  }// End of TEST Area
70 
71 
72 
73  // This is the first row of equation 32.2.4 on the webpage specified.
74  const Numeric constant = ((((j1-j2-m3)%2)==0)?(1.):(-1.)) * sqrt(factorials(
75  MakeArray<Index>( (+j1+j2-j3).toIndex(),
76  (+j1-j2+j3).toIndex(),
77  (-j1+j2+j3).toIndex(),
78  (j1+m1).toIndex(),
79  (j1-m1).toIndex(),
80  (j2+m2).toIndex(),
81  (j2-m2).toIndex(),
82  (j3+m3).toIndex(),
83  (j3-m3).toIndex()),
84  MakeArray<Index>( (J+1).toIndex())));
85 
86  Numeric sum=0; Rational s=0;
87 
88  // This is the sum of the same function.
89  while( (j1+j2-j3-s).toIndex() >= 0 &&
90  (j1-m1-s).toIndex() >= 0 &&
91  (j2+m2-s).toIndex() >= 0 )
92  {
93  if( (j3-j2+m1+s).toIndex() >= 0 &&
94  (j3-j1-m2+s).toIndex() >= 0 )
95  {
96  sum += ((((s)%2)==0)?(1.):(-1.)) * factorials(
97  MakeArray<Index>( 1),
99  (j1+j2-j3-s).toIndex(),
100  (j1-m1-s).toIndex(),
101  (j2+m2-s).toIndex(),
102  (j3-j2+m1+s).toIndex(),
103  (j3-j1-m2+s).toIndex()));
104  }
105  s++;
106  }
107 
108  return constant*sum;
109 }
110 
111 
119 Numeric wigner6j(const Rational j1,const Rational j2,const Rational j3,
120  const Rational J1,const Rational J2,const Rational J3)
121 {
122 
123  if( !( triangular_inequality( j1, j2, j3 ) ) )
124  return 0;
125 
126  if( !( triangular_inequality( j1, J2, J3 ) ) )
127  return 0;
128 
129  if( !( triangular_inequality( J1, j2, J3 ) ) )
130  return 0;
131 
132  if( !( triangular_inequality( J1, J2, j3 ) ) )
133  return 0;
134 
135  //Fancy simplification
136  if(J1.toIndex()==1)
137  if(j2==J3)
138  if(j3 == J2)
139  return 2 * ((((j1+j2+j3+1).toIndex()%2)==0)?(1.):(-1.)) *
140  ( j2.toNumeric()*(j2.toNumeric()+1) + j3.toNumeric()*(j3.toNumeric()+1) - j1.toNumeric()*(j1.toNumeric()+1) ) /
141  sqrt( 2*j2.toNumeric()*(2*j2.toNumeric()+1)*(2*j2.toNumeric()+2)*2*j3.toNumeric()*(2*j3.toNumeric()+1)*(2*j3.toNumeric()+2) );
142 
143  Numeric sum=0; Rational s=0;
144  //std::cout<<std::endl<<"Wigner6j (" << j1 <<" "<<j2 <<" "<<j3<<"; " << J1 <<" " << J2 << " "<<J3<<")"<<std::endl;
145  // Complicated sum
146  while( (j1+j2+J1+J2-s).toIndex() >= 0 &&
147  (j2+j3+J2+J3-s).toIndex() >= 0 &&
148  (j3+j1+J3+J1-s).toIndex() >= 0 )
149  {
150  if( (s-j1-j2-j3).toIndex() >= 0 &&
151  (s-j1-J2-J3).toIndex() >= 0 &&
152  (s-J1-j2-J3).toIndex() >= 0 &&
153  (s-J1-J2-j3).toIndex() >= 0 )
154  {
155  sum += ((((s)%2)==0)?(1.):(-1.))*factorials(
156  MakeArray<Index>( s.toIndex()+1),
157  MakeArray<Index>( (s-j1-j2-j3).toIndex(),
158  (s-j1-J2-J3).toIndex(),
159  (s-J1-j2-J3).toIndex(),
160  (s-J1-J2-j3).toIndex(),
161  (j1+j2+J1+J2-s).toIndex(),
162  (j2+j3+J2+J3-s).toIndex(),
163  (j3+j1+J3+J1-s).toIndex()));
164  }
165  s++;
166  }
167  //std::cout<<"s became " << s-1 << std::endl;
168 
169  const Numeric tc = sqrt(
170  triangle_coefficient(j1,j2,j3) *
171  triangle_coefficient(j1,J2,J3) *
172  triangle_coefficient(J1,j2,J3) *
173  triangle_coefficient(J1,J2,j3)
174  );
175 
176  return tc * sum;
177 }
178 
179 // /*!
180 // * Equation 34.2.5 in http://dlmf.nist.gov/34.2
181 // */
182 Numeric wigner9j(){return NAN;}//Placeholder
183 
184 
186  Rational Jk_lower, Rational Jl_lower,
187  Rational Jk_upper, Rational Jl_upper)
188 {
189  const Numeric A1 = wigner3j(Nl,Nk,L,0,0,0);
190  if( A1 == 0)
191  return 0;
192  //std::cout << "input to wigner3j is: (" << Nl <<", " << Nk << ", " << L << ")\n";
193 
194  const Numeric A2 = wigner6j(L,Jk_upper,Jl_upper,1,Nl,Nk);
195  if( A2 == 0)
196  return 0;
197 
198  const Numeric A3 = wigner6j(L,Jk_lower,Jl_lower,1,Nl,Nk);
199  if( A3 == 0)
200  return 0;
201 
202  const Numeric A4 = wigner6j(L,Jk_upper,Jl_upper,1,Jl_lower,Jk_lower);
203  if( A4 == 0)
204  return 0;
205 
206  //std::cout << "A1="<<A1<<", A2="<<A2 << ", A3="<<A3<<", A4="<<A4<<std::endl;
207 
208  return A1*A2*A3*A4;
209 }
210 
211 
213 {
214  Rational nom1=(a+b-c), nom2=(a-b+c), nom3=(-a+b+c), denom1=(a+b+c+1);
215  return factorials(
216  MakeArray<Index>(nom1.toIndex(),nom2.toIndex(),nom3.toIndex()),
217  MakeArray<Index>(denom1.toIndex()));
218 }
219 
220 bool triangular_inequality(const Rational x, const Rational y, const Rational z)
221 {
222  if( abs(x-y) > z )
223  return false;
224  if( (x+y) < z )
225  return false;
226 
227  return true;
228 }
229 
230 
231 // Helper function that returns all primes lower than the input number.
232 void primes(ArrayOfIndex& output, const Index input)
233 {
234  output.resize(0);
235 
236  if(input<2)
237  return;
238 
239  output.push_back(2);
240 
241  if(input==2) return;
242 
243  Index ii = 3,JJ=1;
244 
245  while( ii<=input )
246  {
247  Index test = 1;
248  for(Index jj=0;jj<JJ;jj++)
249  {
250  // If it is divisible by a previous prime number then it is not a prime number.
251  if((ii%output[jj])==0)
252  {
253  test = 0;
254  break;
255  }
256  }
257  // If it is divisible by a previous prime number then it is not a prime number otherwise it is.
258  if( test )
259  {
260  output.push_back(ii);
261  JJ++;
262  }
263 
264  ii++;
265  }
266 
267 }
268 
269 
270 // Helper function that returns the powers for describing factorial of input
271 void powers(ArrayOfIndex& output, const ArrayOfIndex primes, const Index input)
272 {
273  Index numbers = primes.nelem();
274 
275  output.resize(numbers);
276 
277  for(Index ii=0;ii<numbers;ii++)
278  {
279  output[ii] = 0;
280  Index num = primes[ii];
281  while( input/num != 0 )
282  {
283  output[ii] += input/num;
284  num *= primes[ii];
285  }
286  }
287 }
288 
289 
290 // Function that calculates prod(NomFac!) / prod(DenomFac!)
291 Numeric factorials(const ArrayOfIndex& NomFac, const ArrayOfIndex& DenomFac)
292 {
293  //std::cout<<NomFac<<", "<<DenomFac<<std::endl;
294  Index max_nom_ind, max_denom_ind;
295 
296  {//Test for Errors
297  std::ostringstream os;
298  Index test_if_input_is_crazy = 0;
299 
300  if(NomFac.nelem()>0)
301  {
302  max_nom_ind = max(NomFac);
303  if( min(NomFac)<0 )
304  {
305  test_if_input_is_crazy++;
306  os << "Negative values in the factorial nominator.\n This is: " << NomFac << std::endl;
307  }
308  }
309  else
310  max_nom_ind = 0;
311  if(DenomFac.nelem()>0)
312  {
313  max_denom_ind = max(DenomFac);
314  if( min(DenomFac)<0 )
315  {
316  test_if_input_is_crazy++;
317  os << "Negative values in the factorial denominator.\n This is: " << DenomFac << std::endl;
318  }
319  }
320  else
321  max_denom_ind = 0;
322  if(test_if_input_is_crazy>0)
323  throw std::runtime_error(os.str());
324  }
325 
326  const Index max_ind = (max_denom_ind<max_nom_ind) ? max_nom_ind : max_denom_ind;
327 
328  // All primes are based on the largest input
329  ArrayOfIndex all_primes, all_powers, temp;
330  primes(all_primes, max_ind);
331  all_powers.resize(all_primes.nelem());
332 
333  // Get powers and add to all_powers
334  for( Index jj = 0; jj < NomFac.nelem(); jj++ )
335  {
336  powers(temp,all_primes,NomFac[jj]);
337 
338  for( Index ii = 0; ii < temp.nelem(); ii++ )
339  all_powers[ii] += temp[ii];
340  }
341 
342  // Get powers and subtract from all_powers
343  for( Index jj = 0; jj < DenomFac.nelem(); jj++ )
344  {
345  powers(temp,all_primes,DenomFac[jj]);
346 
347  for( Index ii = 0; ii < temp.nelem(); ii++ )
348  all_powers[ii] -= temp[ii];
349  }
350 
351  Numeric total=1;
352  bool still_possible=true; Index num=all_primes.nelem();
353 
354  //std::cout<<all_primes<<"; "<<all_powers<<std::endl;
355 
356  // This loop might be a time consumer.
357  while( still_possible )
358  {
359  Index test=1;
360  for( Index jj = 0 ; jj < num; jj++)
361  {
362  if(all_powers[jj]>0)
363  {
364  if(total > 1e20 && min(all_powers)<0)
365  continue;
366  total *= (Numeric) all_primes[jj];
367  all_powers[jj]--;
368  }
369  else if(all_powers[jj]<0)
370  {
371  if(total < 1e-20 && max(all_powers)>0)
372  continue;
373  total /= (Numeric) all_primes[jj];
374  all_powers[jj]++;
375  }
376  else
377  {
378  test++;
379  }
380  }
381  if(test>=num)
382  still_possible=false;
383  }
384 
385  return total;
386 }
387 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void powers(ArrayOfIndex &output, const ArrayOfIndex primes, const Index input)
Index nelem() const
Number of elements.
Definition: array.h:176
Explicit construction of Arrays.
Definition: make_array.h:51
Numeric wigner9j()
Numeric wigner6j(const Rational j1, const Rational j2, const Rational j3, const Rational J1, const Rational J2, const Rational J3)
Numeric toNumeric() const
Definition: rational.h:52
#define max(a, b)
Definition: continua.cc:20461
void primes(ArrayOfIndex &output, const Index input)
Index toIndex() const
Definition: rational.cc:32
#define abs(x)
Definition: continua.cc:20458
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
This can be used to make arrays out of anything.
Definition: array.h:40
#define min(a, b)
Definition: continua.cc:20460
Numeric factorials(const ArrayOfIndex &NomFac, const ArrayOfIndex &DenomFac)
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
#define temp
Definition: continua.cc:20773
Numeric ECS_wigner(Rational L, Rational Nl, Rational Nk, Rational Jk_lower, Rational Jl_lower, Rational Jk_upper, Rational Jl_upper)
bool triangular_inequality(const Rational x, const Rational y, const Rational z)
Numeric triangle_coefficient(const Rational a, const Rational b, const Rational c)