function Add_LoadCase(BeamZ_file, new_load_case,BeamZ_file_out,name)
% Add_LoadCase ( BEAMZ_FILE, NEW_LOAD_CASE )
% Add a new load case to the BeamZ data file named BEAMZ_FILE.
% NEW_LOAD_CASE specifies the nodes at which forces and moments are applied
% It should contain the following fields : nodes, fx, fy, mz.
% The lengths of NEW_LOAD_CASE.{nodes, fx, fy, mz} should be the same. Use
% zeros if no load nor moment is applied to some specific nodes.
% This function returns a new data file called as the initial data file,
% plus '_modified'. E.g. if BEAMZ_FILE is 'test.BeamZ',  BEAMZ_FILE_OUT is
% 'test_modified.BeamZ'.
% 
% Add_LoadCase ( BEAMZ_FILE, NEW_LOAD_CASE, BEAMZ_FILE_OUT ) allows to
% specify the new file name. This argument is optional.
%
% Use Add_LoadCase recursively if multiple load cases have to be added (in
% this case, use three arguments and BEAMZ_FILE_OUT shall be the same as
% BEAMZ_FILE, so that it is replaced for every added load case).
%
% Example : in this example, a load case is added in Test_bidon.BeamZ.
% Forces are added at nodes 1 and 3.
% At node 1, forces in the x and y directions are 0 and -1. No moment.
% At node 3, forces in the x and y directions are 1 and 5. Moment = 7.
%
% BeamZ_file = 'Test_bidon.BeamZ';
% new_load_case.name = 'Equiv. Static Load 1';
% new_load_case.nodes = [1 3];
% new_load_case.fx    = [0 1];
% new_load_case.fy    = [-1 5];
% new_load_case.mz    = [0 7];
%

L=load(BeamZ_file,'-mat');


nLoads = length(L.Userdata.Loads);

F_node = [];
M_node = [];
k_f=1;
k_m=1;

for i_node=1:length(new_load_case.nodes)
    
    node = new_load_case.nodes(i_node);

    if node > length(L.Userdata.Nodes.XNOD)
        warning ('Wrong node number !!')
        node = nan;
    end
    
    if new_load_case.fx(i_node) ~=0
        F_node(k_f).node  = node;
        F_node(k_f).loadX = new_load_case.fx(i_node);
        F_node(k_f).loadY = 0;
        F_node(k_f).loadx = NaN;
        F_node(k_f).loady = NaN;
        F_node(k_f).magnitude = new_load_case.fx(i_node);
        F_node(k_f).direction = 'GlobalX';
        k_f=k_f+1;
    end

    if new_load_case.fy(i_node) ~=0
        F_node(k_f).node  = node;
        F_node(k_f).loadX = 0;
        F_node(k_f).loadY = new_load_case.fy(i_node);
        F_node(k_f).loadx = NaN;
        F_node(k_f).loady = NaN;
        F_node(k_f).magnitude = new_load_case.fy(i_node);
        F_node(k_f).direction = 'GlobalY';
        k_f=k_f+1;
    end

    if new_load_case.mz(i_node) ~=0
        M_node(k_m).node  = node;
        M_node(k_m).magnitude = new_load_case.mz(i_node);
        k_m=k_m+1;
    end
    
end

Userdata = L.Userdata;
Userdata.Loads(nLoads+1).F_node = F_node;
Userdata.Loads(nLoads+1).M_node = M_node;


if isempty(name); name='new load case'; end
Userdata.Loads(nLoads+1).Name = name;
Userdata.Loads(nLoads+1).isEditbl = 1;

if nargin<3
    s = strsplit('Test_bidon.BeamZ','.');
    BeamZ_file_out = [s{1} '_modified.BeamZ'];
end

save(BeamZ_file_out,'Userdata');