Расчёт 3j-символов специального вида, когда один из j=1

[madschumacher]

Добрый день.

(Дисклеймер для Админов)

Помещаю тему в ПРР(Ф), т.к. полагаю, что с 3j-символами именно посетители этого раздела имеют чаще всего дело

Прошу прощения за беспокойство у Всех. Вновь у меня идиотические проблемы, связанные с моим непроходимым идиотизмом.
Понадобились мне (на определённый чёрт) Вигнеровские 3j-символы определённого вида:
$ \begin{pmatrix} j & 1 & j' \\ m & \sigma & m' \\ \end{pmatrix} $
Естественно, выражения для ненулевых коэффициентов такого вида доступны в Священной Книге (то бишь 3-м Ландавшице):
$ \begin{pmatrix} j & 1 & j \\ m & 0 & -m \\ \end{pmatrix} = (-1)^{j+m+1} \sqrt{\frac{2m}{2j(2j+1)(2j+2)}} $
$ \begin{pmatrix} j & 1 & j \\ m & 1 & -(m+1) \\ \end{pmatrix} = (-1)^{j+m+1} \sqrt{\frac{2(j-m)(j+m+1)}{(2j+1)(2j+2)(2j+3)}} $
$ \begin{pmatrix} j+1 & 1 & j \\ m & 0 & -m \\ \end{pmatrix} = (-1)^{j-m} \sqrt{\frac{2(j-m+1)(j+m+1)}{(2j+1)(2j+2)(2j+3)}} $
$ \begin{pmatrix} j+1 & 1 & j \\ m & 1 & -(m+1) \\ \end{pmatrix} = (-1)^{j-m} \sqrt{\frac{(j-m)(j-m+1)}{(2j+1)(2j+2)(2j+3)}} $
(конечно, там они даны не совсем в таком виде, я их преобразовал к нужному мне виду, возможно с ошибками).
Ну и конечно, написал C-шную функцию (точнее несколько версий), чтобы их считать. Последняя имеет такой вид:

Используется синтаксис C

double w3j_sym(int j1, int m1, int sigma, int j2, int m2) 
       {double res=0.0;
          
         if(j1==j2)
          {
           if(sigma==0 && m1==-m2 )
            {res=pow(-1.0, j1+m1+1) * 2*m1 / sqrt(2*j1*(2*j1+1)*(2*j1+2));}

           if(sigma== 1 && m2==-(m1+1) )           
            {res=pow(-1.0, j1+m1+1) * sqrt(2*(j1-m1)*(j1+m1+1)) / sqrt(2*j1*(2*j1+1)*(2*j1+2));}

           if(sigma==-1)
            {res=pow(-1.0, j1+j2+1)*w3j_sym(j1, -m1, -sigma, j2, -m2);}

          }


         if(j1==(j2+1))
          {
          if(sigma==0 && m1==-m2)
            {res=pow(-1.0, j2+m1) * sqrt(2*(j2+m1+1)*(j2-m1+1)) / sqrt(2*(2*j2+1)*(2*j2+2)*(2*j2+3));}
           if(sigma== 1 && m2==-(m1+1) )
            {res=pow(-1.0, j2+m1) * sqrt(2*(j2-m1)*(j2-m1+1)) / sqrt(2*(2*j2+1)*(2*j2+2)*(2*j2+3));}
           if(sigma==-1)
            {res=pow(-1.0, j1+j2+1)*w3j_sym(j1, -m1, -sigma, j2, -m2);}
          }

        return res;
       }
 

Косячность реализации (что j1 $\geq$ j2) компенсируется тем, что так они в эту функцию потом и поступают.
Но, сравнивая результаты с сетевым 3j-калькулятором http://www-stone.ch.cam.ac.uk/wigner.shtml, а также результаты уже финального кода на тестовом примере с расчетами по аналогичной проге я получаю грусть, печаль и разочарование в своих умственных способностях.

Прошу помочь кто чем может! (статейкой/книжкой, где расчёты подобных символов разжёваны для таких идиотов, как я, уже готовым кодом или алгоритмом, который можно без угрызений совести вставить в свой, или просто советом и комментарием, что я идиот) Всему буду рад.

Я уже перерыл кучу этих наших Интернетов, как русских, так и заграничных. Посмотрел и другие книги (такие, как Варшалович, Зар и т.д.) и статьи, но от идиота головного мозга ничего не помогает.

Поэтому и надеюсь на помощь живительных советов и пенделей от Здешних Обитателей.

[amon]

Что бы я сделал. Я бы вбил в Математику общее выражение из ЛЛ предварительно сверившись, например, с Варшаловичем на предмет возможных ошибок в ЛЛ, и заставил бы железяку сосчитать мне явные выражения для Вашего случая. После этого, не мудрствуя лукаво, я бы заставил Математику выплюнуть сишный код полученного выражения и copy-paste'ом вставил бы его в свою программу.

[g______d]

В коде в первом варианте $2 m_1$ не под корнем, так и должно быть?

[madschumacher]

g______d, да, действительно ошибся. Должно быть
$\begin{pmatrix} j & 1 & j \\ m & 0 & -m \\ \end{pmatrix} = (-1)^{j+m+1} \frac{2m}{\sqrt{2j(2j+1)(2j+2)}}$

g______d в сообщении #1319193 писал(а):

Я бы вбил в Математику общее выражение из ЛЛ предварительно сверившись, например, с Варшаловичем на предмет возможных ошибок в ЛЛ

К сожалению проблема не сколько в самом вбивании формул, сколько в том, чтобы их размножить, поскольку в этих книгах приводятся минимальные наборы формул, из которых потом при помощи преобразований перестановок и инверсий генерируются все остальные. В первоначальном варианте я все формулы посчитал вручную и забил через if-ы с соответствующими условиями. Но, видимо, неправильно.

А код для C, который генерирует Mathematika 11.1, которая у меня на рабочем компе, для нужных мне 3j-символов, несколько некомпилируемый (он содержит такие штуки, как piecewise и факториалы). Ну ладно, это всё можно закодить, но факториал для $j > 20$ -- это уже жёстко...

[madschumacher]

amon, спасибо за совет. Я взял Wolfram Mathematica 11.1 и посчитал все возможные варианты 3j-символов (правда, пришлось поупрощать ручками, ну да и ладно, вроде не накосячил, я надеюсь). Чтобы не потерять (и если м.б. кому-то понадобится) публикую все эти случаи тут.
Общий вид 3j-символа:
$\begin{pmatrix} j_1 & 1 & j_3 \\ m_1 & m_2 & m_3 \\ \end{pmatrix}$

• $j_1 = j_3 = j$:
  
  • $m_2 = +1$:
    $ \begin{pmatrix} j & 1 & j \\ m_1 & +1 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3+1} \cdot \frac{(-1)^{m_3 - j}}{\sqrt{2}} \cdot \sqrt{ \frac{(1+j+m_3)(j-m_3)}{j (j+1) (2j+1)} } $
  • $m_2 = 0$:
    $ \begin{pmatrix} j & 1 & j \\ m_1 & 0 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3} \cdot {(-1)^{m_3 - j}} \cdot \frac{m_3}{\sqrt{j (j+1) (2j+1)}}$
    
  • $m_2 = -1$:
    $ \begin{pmatrix} j & 1 & j \\ m_1 & -1 & m_3 \\ \end{pmatrix} = \delta_{1, m_1+m_3} \cdot \frac{(-1)^{j-m_3+1}}{\sqrt{2}} \cdot \sqrt{ \frac{(j+m_3)(j+1-m_3)}{j (j+1) (2j+1)} } $
• $j_1 -1 = j_3 = j$:
  
  • $m_2 = +1$:
    $ \begin{pmatrix} j+1 & 1 & j \\ m_1 & +1 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3+1} \cdot \frac{(-1)^{m_3 - j}}{\sqrt{2}} \cdot \sqrt{ \frac{(1+j+m_3)(2+j+m_3)}{ (j+1) (2j+1)(2j+3)} } $
  • $m_2 = 0$:
    $ \begin{pmatrix} j+1 & 1 & j \\ m_1 & 0 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3} \cdot {(-1)^{1- j+m_3}} \cdot \sqrt{\frac{(j+1-m_3)(j+1+m_3)}{(j+1) (2j+1)(2j+3)}}$
    
  • $m_2 = -1$:
    $ \begin{pmatrix} j+1& 1 & j \\ m_1 & -1 & m_3 \\ \end{pmatrix} = \delta_{1, m_1+m_3} \cdot \frac{(-1)^{m_3-j}}{\sqrt{2}} \cdot \sqrt{\frac{(j+1-m_3)(j+2-m_3)}{(j+1) (2j+1)(2j+3)}} $
  
  
• $j_1 = j_3 - 1 = j$:
  
  • $m_2 = +1$:
    $ \begin{pmatrix} j & 1 & j +1\\ m_1 & +1 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3+1} \cdot \frac{(-1)^{j-m_3-1}}{\sqrt{2}} \cdot \sqrt{ \frac{(j-m_3)(1+j-m_3)}{ (j+1) (2j+1)(2j+3)} } $
  • $m_2 = 0$:
    $ \begin{pmatrix} j & 1 & j+1 \\ m_1 & 0 & m_3 \\ \end{pmatrix} = \delta_{0, m_1+m_3} \cdot {(-1)^{ j+1-m_3}} \cdot \sqrt{\frac{(j+1-m_3)(j+1+m_3)}{(j+1) (2j+1)(2j+3)}}$
    
  • $m_2 = -1$:
    $ \begin{pmatrix} j & 1 & j+1 \\ m_1 & -1 & m_3 \\ \end{pmatrix} = \delta_{1, m_1+m_3} \cdot \frac{(-1)^{j+1-m_3}}{\sqrt{2}} \cdot \sqrt{\frac{(j+m_3)(j+1+m_3)}{(j+1) (2j+1)(2j+3)}} $

Также, на в.с. ещё и привожу C-шную функцию для этих формул:

Используется синтаксис C

double w3j_sym(int j1, int m1, int sigma, int j3, int m3) 
       {double res=0.0;
        int j;  

         if(j1==j3 && j1>0)
          {j=j1;

           if(sigma==1 && (m1+m3+1 == 0) )
             {res=pow(-1.0, m3-j)/sqrt(2.0);
              res*=sqrt((1.0+j+m3)*(j-m3)/(j*(1.0+j)*(1.0+2*j)));
             }

           if(sigma==0 && (m1+m3 == 0) )
             {res=pow(-1.0, m3-j)*m3;
              res/=sqrt(j*(j+1.0)*(1.0+2*j));
             }

           if(sigma==-1 && (m1+m3 == 1) )
             {res=pow(-1.0, j-m3+1)/sqrt(2.0);
              res*=sqrt((j+m3)*(j+1.0-m3)/(j*(1.0+j)*(1.0+2*j)));
             }
 
          }


         if(j1==(j3+1))
          {j=j3;

           if(sigma==1 && (m1+m3+1 == 0) )
             {res=pow(-1.0, m3-j)/sqrt(2.0);
              res*=sqrt((1.0+j+m3)*(2.0+j+m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }           

           if(sigma==0 && (m1+m3 == 0) )
             {res=pow(-1.0, 1-j+m3);
              res*=sqrt((1.0+j-m3)*(1.0+j+m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }           

           if(sigma==-1 && (m1+m3 == 1) )
             {res=pow(-1.0, m3-j)/sqrt(2.0);
              res*=sqrt((1.0+j-m3)*(2.0+j-m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }          
          }         


         if(j3==(j1+1))
          {j=j1;

           if(sigma==1 && (m1+m3+1 == 0) )
             {res=pow(-1.0, j-m3-1)/sqrt(2.0);
              res*=sqrt((j-m3)*(1.0+j-m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }           

           if(sigma==0 && (m1+m3 == 0) )
             {res=pow(-1.0, 1+j-m3);
              res*=sqrt((1.0+j-m3)*(1.0+j+m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }           

           if(sigma==-1 && (m1+m3 == 1) )
             {res=pow(-1.0, j+1-m3)/sqrt(2.0);
              res*=sqrt((j+m3)*(1.0+j+m3)/((1.0+j)*(1.0+2*j)*(3.0+2*j)));
             }          
          }         

        return res;
       }

 

[amon]

(Оффтоп)

Математика, если ее правильно использовать и не считать, что она сама статью напишет, вещь полезная, особенно для таких рекордсменов по количеству ошибок на погонный метр формул или текста как я