[gaussian_matrix,gaussian_rot_matrix]=gaussian_frame(center_of_mass_matrix,dihedral)

function [gaussian_matrix,gaussian_rot_matrix]=gaussian_frame(center_of_mass_matrix,dihedral)

N=length(center_of_mass_matrix(:,1,1));

trials=length(center_of_mass_matrix(1,1,:));

for i=1:trials
    
    bc=center_of_mass_matrix(dihedral(3),1:3,i)-center_of_mass_matrix(dihedral(2),1:3,i);
    
    norm_x_axis=bc/norm(bc);
    
    ba=center_of_mass_matrix(dihedral(1),1:3,i)-center_of_mass_matrix(dihedral(2),1:3,i);
    
    y_axis=cross(ba,norm_x_axis);

    norm_y_axis=y_axis/norm(y_axis);
    
    norm_z_axis=cross(norm_x_axis,norm_y_axis);
    
    lambda=[norm_x_axis;norm_y_axis;norm_z_axis];
    
    gaussian_rot_matrix(1:3,1:3,i)=lambda;
    
    gaussian_matrix(1:N,1:3,i)=center_of_mass_matrix(1:N,1:3,i)*lambda';
    
end