Status: Experimental / Incomplete
Origin: src/ISGL-3body/precomputation.jl (former lines 49-51 and 462+)
Reason archived: unfinished spin-orbit precomputation
Dependencies: imaxSO_dict, mijSO_arr_dict, SSO_arr, get_Ds, Si_funSO
Restoration steps: reintroduce call in precompute_ISGL and add tests
Archived on: 2026-05-21

# 49-52:
# necessary only in case of spin-orbit interactions:
if prod(isempty.(so_indices)) == false
    precompute_SSO(SSO_arr,Clmk_arr,Dlmk_arr,JlL_complete,lL_complete,l_complete,imaxSO_dict,mijSO_arr_dict,kmax_dict,gamma_dict,cleb_arr,temp_S,temp_D1,temp_D2)
end


# 462+
## for spin-orbit interaction:
### S_coeff
@views @inbounds function precompute_SSO(SSO_arr,Clmk_arr,Dlmk_arr,JlL_complete,lL_complete,l_complete,imaxSO_dict,mijSO_arr_dict,kmax_dict,gamma_dict,cleb_arr,temp_arr,tempD1,tempD2)
    # returns the OffsetArray S_arr[la,La,lb,Lb,1:maximax]. only 1:imax(la,La,lb,Lb) is filled with nonzeros.
        
    for JlLa in JlL_complete
        for JlLb in JlL_complete
            (abs(JlLa - JlLb) <= 1 <= JlLa + JlLb) == false && continue # only for allowed values of JlLa,JlLb
            for (la,La) in lL_complete
                (abs(la - La) <= JlLa <= la + La) == false && continue # only for allowed values of la,La,JlLa
                for (lb,Lb) in lL_complete
                    (abs(lb - Lb) <= JlLb <= lb + Lb) == false && continue # only for allowed values of lb,Lb,JlLb
                    for v = 1:6
                        imaxSO_dict[((la,La),(lb,Lb),v)] == 0 && continue # skip if no combinations of mijSO exist
                        
                        SSO_arr[la,La,lb,Lb,JlLa,JlLb,v,:] .= Si_funSO(JlLa,JlLb,la,La,lb,Lb,imaxSO_dict[((la,La),(lb,Lb),v)],mijSO_arr_dict[((la,La),(lb,Lb),v)],kmax_dict,gamma_dict,cleb_arr,Clmk_arr,Dlmk_arr,temp_arr,tempD1,tempD2,v)
                    end
                end
            end
        end
    end
    
end


# to assign pairs of Di,Dj based on v
function get_Ds(v, D1, D2, D3, D4)
    if v == 1
        return D1, D2
    elseif v == 2
        return D1, D3
    elseif v == 3
        return D1, D4
    elseif v == 4
        return D2, D3
    elseif v == 5
        return D2, D4
    else  # v == 6
        return D3, D4
    end
end


@views @inbounds function Si_funSO(JlLa,JlLb,la,La,lb,Lb,imax,mij_arr_i,kmax_dict,gamma_dict,cleb_arr,Clmk_arr,Dlmk_arr,temp_arr,D1,D2,v)
    # returns a vector S(la,La,lb,Lb,J)[i=1,i=2,...,i=maximax] only i=1:imax(la,La,lb,LB) is filled with nonzeros
    
    # reset to zero (the same array is used in each function call)
    for i = 1:lastindex(temp_arr)
        temp_arr[i] = 0.0
    end
    
    M_tot = 0 # not really Mtot, but the index? (not rank) of the tensor operator
    MlLa = min(JlLa, JlLb)
    MlLb = MlLa
    
    mamax = min(la,MlLa+La) # untested
    mamin = max(-la,MlLa-La)
    mbmax = min(lb,MlLb+Lb)
    mbmin = max(-lb,MlLb-Lb)
    if mamax < mamin || mbmax < mbmin
        return temp_arr
    end
    
    # MlLa = MlLb+0 is always true

    # for the cases when cleb_arr was not precomputed (mostly, if lmin,Lmin > 0)
    if cleb_arr[1,M_tot,JlLb,MlLb,JlLa] == 0
        clebJba = PartialWaveFunctions.clebschgordan(1,M_tot,JlLb,MlLb,JlLa,MlLa)
    else
        clebJba = cleb_arr[1,M_tot,JlLb,MlLb,JlLa]
    end
    clebjllabJ = (-1)^(-1+JlLa+JlLb)*sqrt(2*JlLa+1) / clebJba # my formula. untested
    
    for ma = mamin:mamax # no easy truncation from Wigner-Eckart
        Ma = MlLa - ma
        abs(Ma) > La && continue
        
        for mb = mbmin:mbmax
            Mb = MlLb - mb
            abs(Mb) > Lb && continue 
            
            cleb_fac = cleb_arr[la,ma,La,Ma,JlLa]*cleb_arr[lb,mb,Lb,Mb,JlLb]*clebjllabJ
            
            for ka=1:kmax_dict[la,ma]
                for Ka=1:kmax_dict[La,Ma]
                    for kb=1:kmax_dict[lb,mb]
                        for Kb=1:kmax_dict[Lb,Mb]
                            # ccleb_fac = RCFAC
                            ccleb_fac = Clmk_arr[la,ma,ka]*Clmk_arr[La,Ma,Ka]*Clmk_arr[lb,mb,kb]*Clmk_arr[Lb,Mb,Kb] * cleb_fac
                            
                            # D_i, the i is only determined by la,ma,ka
                            D1 .= conj.(Dlmk_arr[la,ma,ka,:])
                            D2 .= conj.(Dlmk_arr[La,Ma,Ka,:])
                            D3 = Dlmk_arr[lb,mb,kb,:]
                            D4 = Dlmk_arr[Lb,Mb,Kb,:]
                            
                            D12 = real(transpose(D1)*D2)
                            D13 = real(transpose(D1)*D3)
                            D14 = real(transpose(D1)*D4)
                            D23 = real(transpose(D2)*D3)
                            D24 = real(transpose(D2)*D4)
                            D34 = real(transpose(D3)*D4)
                            
                            # cleaner call, moved if-tree to extra function get_Ds
                            Di, Dj = get_Ds(v, D1, D2, D3, D4)                            
                            Dcrossz = im*(Di[1]*Dj[2]-Di[2]*Dj[1]) # DixDjy - DiyDjx
                            
                            # calculation of S for all values of i, so it can be accessed later in the sum over i in fillTVS.
                            for i = 1:imax # loop over i inside, as the rest can be done independently
                                m12 = mij_arr_i[i][1] # automatically has the correct v
                                m13 = mij_arr_i[i][2]
                                m14 = mij_arr_i[i][3]
                                m23 = mij_arr_i[i][4]
                                m24 = mij_arr_i[i][5]
                                m34 = mij_arr_i[i][6]
                                
                                # ref: mij_fac = factno:
                                mij_fac = 1/(gamma_dict[m12+1.0]*gamma_dict[m13+1.0]*gamma_dict[m14+1.0]*gamma_dict[m23+1.0]*gamma_dict[m24+1.0]*gamma_dict[m34+1.0])
                                
                                sss = D12^m12*D13^m13*D14^m14*D23^m23*D24^m24*D34^m34*ccleb_fac*mij_fac*Dcrossz
                                temp_arr[i] += sss
                            end
                            
                        end
                    end
                end
            end
        end
    end
    
    return temp_arr
end