
% creates reduced transition density matrix elements from signed connections


function []=trDens(connMbState,aPlusaMinus,sign,sBasen,toRedIndex,redSBasen,Jf,Ji,ci,cf,antalPartiklar,Nmaximal,paritet,Msumma)



Jmin=abs(Jf-Ji);
Jmax=abs(Jf+Ji);
JlistSUPER=[Jmin:1:Jmax];
save={};
for J=JlistSUPER
    save{J+1}=zeros(length(redSBasen),length(redSBasen));
end

minJ=min([Jf Ji]);
Mi=minJ;
Mf=Mi;



for i=1:length(connMbState)
    
    % get amplitudes 
    koppling=ci(connMbState(i,1))*cf(connMbState(i,2))*sign(i);
    
    % get single-particle states
    aPlus=aPlusaMinus(i,1);
    aMinus=aPlusaMinus(i,2);
    
    m=sBasen(aMinus,4);
    j1=sBasen(aMinus,3);
    
    j2=sBasen(aPlus,3);
    
    % gets index from reduced basis
    aRedMinus=toRedIndex(aMinus);
    aRedPlus=toRedIndex(aPlus);
    
    % Compute the relevant J:s
    Jmin=abs(j1-j2);
    Jmax=abs(j1+j2);
    Jlist=[Jmin:1:Jmax];
    Jlist = intersect(Jlist,JlistSUPER);
    
    M=0;
    
    
    for J=Jlist
        temp=save{J+1};
        
        % vi summerar upp detta matriselement

        temp(aRedPlus,aRedMinus)=temp(aRedPlus,aRedMinus)+(-1)^(j2-m)*ClebschGordan(j1,j2,J,m,-m,M) *koppling;
        
        % d?r tecknet p? matris(i,j)=<my|a+a-|ny> ber?knats p? FPGAn
        save{J+1}=temp;
    end
    
    
    
end


% multiplicate with common factor
for J=JlistSUPER
    save{J+1}=sqrt(2*J+1)*save{J+1}./((-1)^(Ji-Mi)*ClebschGordan(Jf,Ji,J,Mf,-Mi,M));
    % alla element i matrisen save{J}
end


% save to file
str=strcat('./Results/output_A=',num2str(antalPartiklar),'Nmax=', ...
    num2str(Nmaximal),'Par=',num2str(paritet),'M=',num2str(Msumma));
fid = fopen(strcat(str,'.txt'),'w');
    
fprintf(fid,str);
fprintf(fid,'\n'); 
fprintf(fid,'\n'); 
for i=1:length(redSBasen)
    fprintf(fid,'#%g',i);
    fprintf(fid,'\t');
    fprintf(fid,'N=%g',redSBasen(i,1));
    fprintf(fid,'\t');
    fprintf(fid,'l=%g',redSBasen(i,2));
    fprintf(fid,'\t');
    fprintf(fid,'j=%g',redSBasen(i,3));
    fprintf(fid,'\n');  
end

fprintf(fid,'\n');
fprintf(fid,'\n');


for J=Jlist
    fprintf(fid,'Jtrans%d',J);
    fprintf(fid,'\n');
    matrix=save{J+1};
    
    for i=1:length(matrix)

        for j=1:length(matrix)
            if matrix(i,j)~=0 % Hittat output
                fprintf(fid,'a+=%g',i);
                fprintf(fid,'\t');
                fprintf(fid,'a-=%g',j);
                fprintf(fid,'\t');
                fprintf(fid,'trDens=%g',matrix(i,j));
                fprintf(fid,'\n');
            end
        end
    end

    fprintf(fid,'\n');
end
fclose(fid);


end