clc; clear; close all

% PARAMETRIC STUDY FOR TWIN-CARCASS WHEEL
% GEOMETRIC PARAMETERS

g = 1.62;

R_ext_vec= 0.5;                                      % radius
B_vec= 0.4;                                          % width
t_ext_vec= 0.003;                                    % outer ring thickness
t_int_vec= 0.003;                                    % inner ring thickness
t_rinf_vec= 0.003;                                   % reinforcement thickness 
N_vec= 8;                                            % number of wheels
M_vec= 7000;                                         % vehicle mass

%%% Variable parameters for the tolerance robustness analysis
% t_rinf_vec= [0.0027,0.00285,0.003,0.00315,0.0033];
% t_int_vec= [0.0027,0.00285,0.003,0.00315,0.0033];
% t_ext_vec= [0.0027,0.00285,0.003,0.00315,0.0033];

% ARCHITECTURAL PARAMETERS 

R_int_vec= R_ext_vec - B_vec./2;
R_ext_cil_vec= R_ext_vec;
N_elem_vec= 9;                         % number of reinforcement
% floor(2*pi*R_int_vec/(R_ext_vec/10));

% MATERIAL PARAMETERS
% Al
% E = 72400E6;
% ni=0.3;
% rho=2700;
% sigma_yAl=415*10^6; % [Pa]

% AlSi10Mg
E=71000E6;
ni=0.33;
rho=2680;
sigma_yAl=264E6;

% MESH PARAMETERS

mesh_vec=R_int_vec/10;

% MATERIAL STRUCTURE

% Aluminium
param_mat.Al.name='Aluminium';
param_mat.Al.prop= [E,ni,rho,sigma_yAl];

% Steel
param_mat.Ac.name= 'Acciaio';
param_mat.Ac.prop=[E,ni,rho,sigma_yAl];

% Metallic alloy
param_mat.Pol.name= 'AlSi10Mg';
param_mat.Pol.prop=[E,ni,rho,sigma_yAl];

% sigma_y= max(sigma_yAl,sigma_yAl,sigma_yAl);
sigma_y=sigma_yAl;

mat= fieldnames(param_mat); % create a list of materials
n_mat= numel(mat); % count the number of materials

for j = 1:n_mat
    list_mat(j) = param_mat.(mat{j});  %#ok<SAGROW> % create a table with materials and properties
end

%%% Parametric loop

tic;

n_righe = length(R_ext_vec) * length(B_vec) * length(N_vec) * length(t_ext_vec) ...
    *length(t_int_vec)*length(t_rinf_vec)*length(N_elem_vec)*length(M_vec)*length(N_vec); % number of rows in the output matrix
n_colonne= 15;

MATRIX= zeros(n_righe,n_colonne);

totalIter= n_righe;
wait_bar = waitbar(0, 'Starting analysis...');

count=1;
doneIter = 0;

for r= 1:length(R_ext_vec)
    for b= 1:length(B_vec)
        for v= 1: length (t_ext_vec)
            for l=1:length(t_int_vec)
                for g=1:length(t_rinf_vec)
                    for s=1:length(N_elem_vec)
                        for m=1:length(M_vec)
                            for n=1:length(N_vec)
                                % for h=1:length(R_int_vec)

                                R_ext= R_ext_vec(r);
                                R_int=R_int_vec(r);
                                R_ext_cil=R_ext_cil_vec(r);
                                B=B_vec(b);
                                t_ext= t_ext_vec(v);
                                t_int=t_int_vec(l);
                                t_rinf =t_rinf_vec(g);
                                N_elem=N_elem_vec(s);
                                M= M_vec(m);
                                N= N_vec(n);

                                mesh= mesh_vec(r);

                                N_elem_tot= floor(2*pi*R_int_vec/(R_int_vec/5))-1;

                                % SESSION FILE CREATION

                                file_ses= "architecture_twin_carcasses_cilinder.ses";
                                file_sesw=fopen(file_ses,"w");
                                f_db_jou= "architecture_twin_carcasses_cilinder.db.jou";
                                f_conf_ini= "patran.conf.ini";
                                
                                f_ses01='patran.ses.01';
                                if exist(f_ses01, 'file')
                                    delete(f_ses01);
                                end
                                
                                if exist(f_db_jou, 'file')
                                    delete(f_db_jou);
                                end
                                
                                if exist(f_conf_ini, 'file')
                                    delete(f_conf_ini);
                                end
                                
                                fullpath = mfilename('fullpath');
                                path_App = strrep(fileparts(fullpath), '\', '/');
                                path_Pat = "C:\Program Files\MSC.Software\NaPa_SE\20241\Patran\bin\patran.exe";
                                model_path = strrep(fullfile(path_App, "architecture_twin_carcasses_cilinder.db"), '\', '/');   
                                [path1, ~, ~] = fileparts(path_Pat);
                                [path_folder_pat, ~, ~] = fileparts(path1);
                                path_folder_pat = strrep(path_folder_pat, '\', '/');
                                template_path = strrep(fullfile(path_folder_pat, 'template.db'), '\', '/');
                                
                                fprintf(file_sesw, 'uil_file_new.go( "%s", "%s" )\n', template_path, model_path);
                                fprintf(file_sesw, 'set_current_dir( "%s" )\n', path_App );


                                % ARCHITECTURE STRUCTURE

                                W= M*g; % total weight [N]
                                P= W/N; % Vertical force on a wheel
                                
                                param_str= struct ('R_ext_cil',R_ext_cil, ...
                                    'R_int',R_int, ...
                                    'B', B, ...
                                    't_ext', t_ext, ...
                                    't_int',t_int,...
                                    't_rinf',t_rinf,...
                                    'W',W, ...
                                    'P', P, ...
                                    'N_elem',N_elem, ...
                                    'N_elem_tot',N_elem_tot,...
                                    'mesh',mesh);
                                
                                
                                % PROPERTY STRUCTURE
                                
                                % Shell
                                param_Shell.name = "Shell_" + string(mat)';
                                param_Shell.th.name = ["t_int", "t_rinf", "t_ext"];
                                param_Shell.th.val= [t_int,t_rinf,t_ext];

                                % SESSION FILE PARAMETERS

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ GEOMETRY PARAMETERS\n');
                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, 'REAL R_ext_cil = %d\n', param_str.R_ext_cil);
                                fprintf(file_sesw, 'REAL R_int = %d\n', param_str.R_int);
                                fprintf(file_sesw, 'REAL B = %d\n', param_str.B);
                                fprintf(file_sesw, 'REAL t_ext = %d\n', param_str.t_ext);
                                fprintf(file_sesw, 'REAL t_int = %d\n', param_str.t_int);
                                fprintf(file_sesw, 'REAL t_rinf = %d\n', param_str.t_rinf);
                                fprintf(file_sesw, 'REAL P = %d\n', param_str.P);
                                fprintf(file_sesw, 'INTEGER N_elem = %i\n', param_str.N_elem);
                                
                                fprintf(file_sesw, 'REAL mesh = %d\n', param_str.mesh);

                                % SESSION FILE GEOMETRY

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ GEOMETRY\n');
                                fprintf(file_sesw, newline);

                                % hub
                                
                                p_id=0;
                                c_id=0;
                                s_id=0;

                                % n_an= mod((0:360/param_str.N_elem:360)+270,360);  % point indices 
                                n_an = round(mod(linspace(270, 270+360, param_str.N_elem+1), 360));
             
                                for i=1:param_str.N_elem

                                    fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n'); 
                                    fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, %d, %d., "Coord 0.3", "[0 0 0]", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_int,n_an(i)-360/param_str.N_elem_tot/2,n_an(i));
                                    c_id=c_id+1;
                                    fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, %d, %d., "Coord 0.3", "[0 0 0]", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_int,n_an(i),n_an(i)+360/param_str.N_elem_tot/2);   
                                    c_id=c_id+1;

                                    p_id=p_id+3;

                                end
                                
                                list_p_tot=1:p_id;
                                list_p_mean= 2:3:list_p_tot(end);
                                p_id= p_id-length(list_p_mean);
                                list_p1=setdiff(list_p_tot,list_p_mean);
                                list_p2=setdiff(list_p_tot,list_p_mean); % the one I need

                                for i=1:param_str.N_elem % it also creates n_elem useless points

                                    if i==param_str.N_elem
                                        fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'sgm_const_curve_2d_arc2point_v2( "%i", 2, %d, TRUE, FALSE, 1, "Coord 0.3", "", "Point %i", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_int,list_p1(end),list_p2(1) );
                                        c_id=c_id+1;
                                        p_id=p_id+1;

                                    elseif i<param_str.N_elem

                                        fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'sgm_const_curve_2d_arc2point_v2( "%i", 2, %d, TRUE, FALSE, 1, "Coord 0.3", "", "Point %i", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_int,list_p1(i+1),list_p1(i+2) );
                                        c_id=c_id+1;
                                        p_id=p_id+1;
                                        list_p1= setdiff(list_p1,list_p1(1));

                                    end

                                end
                                
                                p_id_right=p_id+param_str.N_elem;
                                list_p3= [setdiff(list_p_tot,list_p_mean),p_id+1:p_id_right];
                                fprintf(file_sesw, 'STRING asm_delete_point_deleted_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'asm_delete_point( "Point %i:%i", asm_delete_point_deleted_ids )\n',list_p3(end-param_str.N_elem+1),list_p3(end));
                                % p_id=p_id-param_str.N_elem;
                                p_id_right=p_id_right-param_str.N_elem;

                                list_c=1:c_id; % the one I need

                                fprintf(file_sesw, 'STRING sgm_sweep_surface_e_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_surface_extrude( "%i", "<0 0 %d>", 1., 0., "[0 0 0]", "Coord 0", "Curve %i:%i", sgm_sweep_surface_e_created_ids )\n',s_id+1,param_str.B/2,1,c_id);
                                s_id=s_id+param_str.N_elem*2+param_str.N_elem;
                                % p_id=p_id+param_str.N_elem*2+param_str.N_elem;
                                p_id_right=p_id_right+param_str.N_elem*2+param_str.N_elem;
                                list_p4=list_p2(end)+1:p_id_right;

                                fprintf(file_sesw, 'sgm_const_surface_extrude( "%i", "<0 0 %d>", 1., 0., "[0 0 0]", "Coord 0", "Curve %i:%i", sgm_sweep_surface_e_created_ids )\n',s_id+1,param_str.B/2+param_str.B/2,1,c_id);
                                s_id=s_id+param_str.N_elem*2+param_str.N_elem;
                                % p_id=p_id+param_str.N_elem*2+param_str.N_elem;
                                p_id_right=p_id_right+param_str.N_elem*2+param_str.N_elem;
                                list_p5=list_p4(end)+1:p_id_right;

                                list_s= 1:s_id;
                                list_p= 1:p_id_right;

                                % REINFORCEMENTS

                                % create reinforcement curves

                                fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_curve_2d_arc2point_v2( "%i", 2, %d, TRUE, FALSE, 2, "Coord 0.1", "", "Point %i", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.B/2/2,list_p_mean(1),list_p4(2)); % B*3/8 /2
                                fprintf(file_sesw, 'sgm_const_curve_2d_arc2point_v2( "%i", 2, %d, TRUE, FALSE, 2, "Coord 0.1", "", "Point %i", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+2,param_str.B/2/2,list_p4(2),list_p5(2)); % B*3/8 /2          
                                c_id=c_id+2;
                                p_id_right= p_id_right+2;

                                list_c_rinf1=list_c(end)+1:c_id;
                                list_p_rinf1=list_p(end)+1:p_id_right; % useless

                                % create surface from curves

                                fprintf(file_sesw, 'STRING sgm_sweep_surface_r_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_surface_revolve( "%i", "Coord 0.3", %i, 0, "Coord 0", "Curve %i %i", sgm_sweep_surface_r_created_ids )\n',s_id+1,360/param_str.N_elem_tot/2,c_id-1,c_id);
                                s_id=s_id+2;
                                fprintf(file_sesw, 'sgm_const_surface_revolve( "%i", "Coord 0.3", %i, 0, "Coord 0", "Curve %i %i", sgm_sweep_surface_r_created_ids )\n',s_id+1,-360/param_str.N_elem_tot/2,c_id-1,c_id);
                                s_id=s_id+2;

                                list_s_rinf1= list_s(end)+1:s_id;
                              
                                % list_s_rinf_tot= list_s(end)+1:s_id;                             
                                
                                % create outer diameter

                                fprintf(file_sesw, 'STRING asm_grid_extract_cu_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'asm_const_grid_extract_v1( "%i", "Curve %i %i", 0.5, 1, asm_grid_extract_cu_created_ids )\n',p_id_right+1,c_id-1,c_id);
                                p_id_right=p_id_right+2;

                                fprintf(file_sesw, 'STRING sgm_sweep_curve_rev_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_curve_revolve( "%i", "Coord 0.3", %d, 0., "Coord 0", "Point %i %i", sgm_sweep_curve_rev_created_ids )\n',c_id+1,360/param_str.N_elem,p_id_right-1,p_id_right);
                                c_id=c_id+2;
                                p_id_right=p_id_right+2;

                                fprintf(file_sesw, 'STRING sgm_sweep_surface_e_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_surface_extrude( "%i", "<0 0 %d>", 1., 0., "[0 0 0]", "Coord 0", "Curve %i %i", sgm_sweep_surface_e_created_ids )\n',s_id+1,param_str.B/4,c_id-1,c_id);
                                s_id=s_id+2;
                                p_id_right=p_id_right+4;

                                fprintf(file_sesw, 'STRING sgm_sweep_surface_e_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_const_surface_extrude( "%i", "<0 0 %d>", 1., 0., "[0 0 0]", "Coord 0", "Curve %i %i", sgm_sweep_surface_e_created_ids )\n',s_id+1,-param_str.B/4,c_id-1,c_id);
                                s_id=s_id+2;
                                p_id_right=p_id_right+2;
                                list_p6=list_p5(end)+1:p_id_right; 

                                list_s_rinf2= list_s_rinf1(end)+1:s_id;

                                % repeat the surfaces

                                fprintf(file_sesw, 'STRING sgm_transform_surf__created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'sgm_transform_rotate( "", "surface", "Coord 0.3", %i, 0., "Coord 0", %i, FALSE, "Surface %i:%i", sgm_transform_surf__created_ids )\n', 360/param_str.N_elem, param_str.N_elem-1, list_s_rinf1(1),s_id);
                                s_id=s_id+8*(param_str.N_elem-1);
                                p_id_right=p_id_right+4*(param_str.N_elem-1)-1;

                                list_s_rinf_tot= list_s_rinf1(1):s_id;

                                % create new cylindrical reference system

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, 'STRING asm_create_cord_3po_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'asm_const_coord_3point( "1", "Coord 0", 2, "[0 0 0]", "[0 0 1]", "[1 0 0]", asm_create_cord_3po_created_ids )\n');                              

                                % Create force application point (at the wheel center)

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, 'STRING asm_create_grid_xyz_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'asm_const_grid_xyz( "%i", "[ 0 0 %d]", "Coord 1", asm_create_grid_xyz_created_ids )\n',p_id_right+1,param_str.B/2);
                                p_id_right=p_id_right+1;

                                % MATERIAL

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ MATERIAL\n');
                                fprintf(file_sesw, newline);

                                for j = 1:n_mat
                                    name_mat = list_mat(j).name;
                                    prop_mat = list_mat(j).prop;
                                    fprintf(file_sesw, ['material.create( "Analysis code ID", 1, "Analysis type ID", 1, "%s", 0, "Date: 21-Nov-25 ' ...
                                        '          Time: 16:49:58", "Isotropic", 1, "Directionality", 1, "Linearity", 1, "Homogeneous", 0, "Linear Elastic", 1,' ...
                                        ' "Model Options & IDs", ["", "", "", "", ""], [0, 0, 0, 0, 0], "Active Flag", 1, "Create", 10, "External Flag", FALSE, ' ...
                                        '"Property IDs", ["Elastic Modulus", "Poisson Ratio", "Density"], [2, 5, 16, 0], "Property Values", ["%d", "%d", "%d", "0"] )\n'], ...
                                        name_mat, prop_mat(1), prop_mat(2), prop_mat(3));
                                end

                                % PROPERTIES

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ PROPERTY\n');
                                fprintf(file_sesw, newline);

                                %%% case t_int=Al, t_rinf=Ac 

                                selected_mat= [3 3]; % 3];

                                name_th = param_Shell.th.name;
                                n_th = length(param_Shell.th.val);  % there are 3 different thicknesses
                                n_sup = {list_s, list_s_rinf_tot}; %, list_s_rinf};   % surface lists

                                for k=1:n_th-1 % there are 2 different thicknesses in this case

                                    name_case_shell=  sprintf("Shell_case1_%s", name_th{k});
                                    material_name = list_mat( selected_mat(k) ).name;
                                    val_th = param_Shell.th.val(k);
                                    list_sup  = n_sup{k};

                                    fprintf(file_sesw, ['elementprops_create( "%s", 51, 25, 35, 1, 1, 20,' ...
                                        ' [13, 20, 36, 4037, 4111, 4118, 4119, 8111, 4401, 4402, 4403, 4404, 4405,' ...
                                        ' 4406, 4407, 4408, 4409], [5, 9, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1, 1, 1, 1, 4, 4],' ...
                                        ' ["m:%s", "", "%d", "", "", "", "", "", "", "", "", "", "", "", "", "", ""],' ...
                                        ' "Surface %i:%i" )\n'], ...
                                        name_case_shell, material_name, val_th, list_sup(1), list_sup(end));

                                end 

                                % LOADS (I apply a point force at the center of the wheel that simulates the
                                % connection with the vehicle and distribute it through RBE3 MPCs to
                                % all mesh nodes on the inner hub surfaces)

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ LOADS\n');
                                fprintf(file_sesw, newline);

                                % create curve for clamping

                                fprintf(file_sesw, 'STRING asm_line_2point_created_ids[VIRTUAL]\n');
                                fprintf(file_sesw, 'asm_const_line_2point( "%i", "Point %i", "Point %i", 0, "", 50., 1, asm_line_2point_created_ids )\n', c_id+1,list_p6(end-3),list_p6(end-1));
                                c_id=c_id+1;

                                fprintf(file_sesw, 'loadsbcs_create2( "Incastro", "Displacement", "Nodal", "", "Static", ["Curve %i"], "Geometry", "Coord 0", "1.", ["< 0 0 0 >", "<  0 0 0   >", "<     >", "<     >"], ["", "", "", ""] )\n', c_id);
                                fprintf(file_sesw, 'loadsbcs_create2( "Force_g", "Force", "Nodal", "", "Static", ["Point %i"], "Geometry", "Coord 0", "1.", ["<    0, %d, 0 >", "<     >", "<     >", "<     >"], ["", "", "", ""] )\n', p_id_right,-param_str.P);


                                % MESH & EQUIVALENCE

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ MESH SURFACE\n');
                                fprintf(file_sesw, newline);

                                fprintf(file_sesw, 'INTEGER fem_create_mesh_surfa_num_nodes\n');
                                fprintf(file_sesw, 'INTEGER fem_create_mesh_surfa_num_elems\n');
                                fprintf(file_sesw, 'STRING fem_create_mesh_s_nodes_created[VIRTUAL]\n');
                                fprintf(file_sesw, 'STRING fem_create_mesh_s_elems_created[VIRTUAL]\n');
                                fprintf(file_sesw, 'fem_create_mesh_surf_4( "IsoMesh", 49152, "Surface 1:%d", 1, ["%d"], "Quad4", "#", "#", "Coord 0", "Coord 0", fem_create_mesh_surfa_num_nodes, fem_create_mesh_surfa_num_elems, fem_create_mesh_s_nodes_created, fem_create_mesh_s_elems_created )\n', list_sup(end), mesh);

                                fprintf(file_sesw, 'REAL fem_equiv_all_x_equivtol_ab\n');
                                fprintf(file_sesw, 'INTEGER fem_equiv_all_x_segment\n');
                                fprintf(file_sesw, 'fem_equiv_all_group5( [" "], 0, "", 1, 1, 0.0049999999, FALSE, fem_equiv_all_x_equivtol_ab, fem_equiv_all_x_segment, TRUE )\n');

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ MESH NODE\n');
                                fprintf(file_sesw, newline);

                                % create dependent node (slave)
                                fprintf(file_sesw, 'STRING fem_create_nodes__nodes_created[VIRTUAL]\n');
                                fprintf(file_sesw, 'fem_create_nodes_1( "Coord 0", "Coord 0", 3, "%d", "Point %i", fem_create_nodes__nodes_created )\n',1000000,p_id_right);

                                % create list A to define all independent nodes (master)
                                fprintf(file_sesw, 'STRING uil_list_create_current_list[VIRTUAL]\n');
                                fprintf(file_sesw, 'list_create_node_ass_geo( "Surface %i:%i", "lista", uil_list_create_current_list )\n',list_s(1),list_s(end));

                                % create RBE3 MPC
                                fprintf(file_sesw, 'fem_create_mpc_nodal2( 1, "RBE3", 0., 2, [TRUE, FALSE], ["0.", "1.0"], ["Node %i", "`lista`"], ["UX,UY,UZ,RX,RY,RZ", "UX,UY,UZ,RX,RY,RZ"] )\n',1000000);

                                % ANALYSIS DECK 

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ ANALYSIS\n');
                                fprintf(file_sesw, newline);


                                fprintf(file_sesw, 'jobfile.open( "architecture_twin_carcasses_cilinder", "ANALYZE NO JOBFILE" )\n');
                                fprintf(file_sesw, 'msc_delete_old_files( "architecture_twin_carcasses_cilinder", ".bdf", ".h5" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "TRANSLATOR", "pat3nas" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "DATABASE", "C:\\Users\\Sofia\\Desktop\\TESI MAGISTRALE\\Matlab tesi\\Architetture ruota\\architecture_twin_carcasses_cilinder\\architecture_twin_carcasses_cilinder.db" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "JOBNAME", "architecture_twin_carcasses_cilinder" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "ANALYSIS TITLE", "MSC.Nastran job created on 23-Nov-25 at" // " 16:49:30" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "ANALYSIS SUBTITLE", "" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "ANALYSIS LABEL", "" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "", "" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "OBJECT", "Entire Model" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "METHOD", "Analysis Deck" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "", "" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "MODEL SUFFIX", ".bdf" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "RESULTS SUFFIX", ".h5" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "SOLUTION TYPE", "LINEAR STATIC" )\n');
                                fprintf(file_sesw, 'jobfile.writei( "SOLUTION SEQUENCE", 101 )\n');
                                fprintf(file_sesw, 'jobfile.writec( "DATABASE RUN", "ON" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "DATA OUTPUT", "PRINT+HDF5" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "OUTPUT2 REQUESTS", "P3 Built In" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "OUTPUT2 FORMAT", "Binary" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "HDF5OUT PRECISION", "HDF5P 32-BIT" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "DIVISION TOLERANCE", "1.0e-08" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "NUMERICAL TOLERANCE", "1.0e-04" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "WRITING TOLERANCE", "1.0e-21" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "GEOM CHECK", "INFORM" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "SHELL NORMAL TOLERANCE", "0.1" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "MASS CALCULATION", "Lumped" )\n');

                                fprintf(file_sesw, 'jobfile.writec( "BULK WRITE", "ON" )\n' );
                                fprintf(file_sesw, 'jobfile.writei( "BULK INPUT 0", 1 )\n' );
                                fprintf(file_sesw, 'jobfile.writec( "BULK INPUT 1", "PARAM,GRDPNT,0" )\n' );
                                fprintf(file_sesw, 'jobfile.writec(  "", "END" )\n' );


                                fprintf(file_sesw, 'jobfile.close()\n');

                                fprintf(file_sesw, 'mscnastran_job.associate_subcases( "101", "architecture_twin_carcasses_cilinder", 1, ["Default"] )\n');
                                fprintf(file_sesw, 'analysis_submit_2( "MSC.Nastran", "architecture_twin_carcasses_cilinder" )\n');

                                str = sprintf('output = system(''"%s" -sfp %s -b -ans yes'');', path_Pat, file_ses);
                                eval(str);


                                % NASTRAN ANALYSIS --> obtain h5 file --> extract maximum
                                % deformation and stress (von Mises) data.

                                path_Nas = "C:\Program Files\MSC.Software\NaPa_SE\20241\Nastran\bin\nastran.exe";
                                file_bdf = "architecture_twin_carcasses_cilinder.bdf";

                                str = sprintf('output=system(''"%s" scr=yes %s'');', path_Nas, file_bdf);
                                eval(str);


                                % Find displacement in y (deltaY) and
                                % wheel stiffness

                                file_h5= 'architecture_twin_carcasses_cilinder.h5';

                                path_disp = '/NASTRAN/RESULT/NODAL/DISPLACEMENT';
                                disp_d = h5read(file_h5,path_disp);
                                U_Y= max(abs(disp_d.Y));

                                K= P/U_Y;

                                % Find Von Mises stress 

                                path_stress = '/NASTRAN/RESULT/ELEMENTAL/STRESS/QUAD_CN';
                                disp_s = h5read(file_h5,path_stress);

                                sigmaxx= disp_s.X1(1,:);
                                sigmayy= disp_s.Y1(1,:);
                                tauxy= disp_s.TXY1(1,:);
                                S_Mises1= max(sqrt(sigmaxx.^2+sigmayy.^2-sigmayy.*sigmaxx+3.*tauxy.^2));

                                sigmaxx2= disp_s.X2(1,:);
                                sigmayy2= disp_s.Y2(1,:);
                                tauxy2= disp_s.TXY2(1,:);
                                S_Mises2= max(sqrt(sigmaxx2.^2+sigmayy2.^2-sigmaxx2.*sigmayy2+3.*tauxy2.^2));

                                S_Mises_max= max(S_Mises1,S_Mises2);

                                DI = 1-S_Mises_max/sigma_y; % Durability index 

                                % Find mass

                                path_mass = '/NASTRAN/RESULT/NODAL/GRID_WEIGHT';
                                disp_m = h5read(file_h5,path_mass);
                                mass= disp_m.MX;


                                %%% RESULTS

                                MATRIX(count, 1) = param_str.R_ext_cil; 
                                MATRIX(count, 2) = param_str.R_int;
                                MATRIX(count, 3) = param_str.B;
                                MATRIX(count, 4) = param_str.t_ext;
                                MATRIX(count, 5) = param_str.t_int;
                                MATRIX(count, 6) = param_str.t_rinf;
                                MATRIX(count, 7) = param_str.W;
                                MATRIX(count, 8) = param_str.P;
                                MATRIX(count, 9) = param_str.N_elem;
                                MATRIX(count, 10) = selected_mat(1); 
                                MATRIX(count, 11) = selected_mat(2);
                                MATRIX(count, 12) = U_Y;
                                MATRIX(count, 13) = S_Mises_max;
                                MATRIX(count, 14) = mass;
                                MATRIX(count, 15) = K;
                                MATRIX(count, 16) = DI;

                                count=count+1;
                                doneIter = doneIter + 1;

                                % Progress bar update

                                elapsed = toc;
                                timePerIter = elapsed / doneIter;
                                remaining = (totalIter - doneIter) * timePerIter;

                                msg = sprintf('Completamento: %.1f \\%% | Tempo rimanente: %.1f s', ...
                                    100*doneIter/totalIter, remaining);

                                waitbar(doneIter/totalIter, wait_bar, msg);

                            end
                        end
                    end
                end
            end
        end
    end
end

close(wait_bar);   % Close the waitbar


MATRIX(count:end,:)= [];

colNames= {'R ext [m]','R int [m]','Width [m]','t ext [m]', 't int [m]', 't rinf [m]',...
    'Vehicle weight [N]','Load [N]','N° rinf','Material ext & int', ...
    'Material rinf','Vertical deflection [m]','Maximum stress [Pa]',...
    'Wheel mass [kg]','Wheel stiffness [kg/s^2]', 'Durability index'};

TAB = array2table(MATRIX, 'VariableNames', colNames); 

%%% SAVE TAB FOR  ARCHITECTURE COMPARISON ANALYSIS
% folderName = 'C:\Users\Sofia\OneDrive - Politecnico di Torino\TESI MAGISTRALE\Matlab tesi\Cartella risultati analisi';
% TAB_all_twin_new= TAB;
% fileName = fullfile(folderName, 'TAB_all_twin_new.mat');
% 
% save(fileName, 'TAB_all_twin_new')


%%% SAVE TAB FOR THE TOLERANCE ROBUSTNESS
% folderName = 'C:\Users\Sofia\OneDrive - Politecnico di Torino\TESI MAGISTRALE\Matlab tesi\Cartella risultati analisi\Tolerance_robustness';
% 
% fileName = fullfile(folderName, 'TAB_tol_twin.mat');
% 
% % 1) If it exists, load the old table
% if isfile(fileName)
%     S = load(fileName, 'TAB_tol_twin');
%     TAB_tol_twin = S.TAB_tol_twin;
% 
%     % 2) Append the new simulation
%     TAB_tol_twin = [TAB_tol_twin; TAB];
% 
%     % 3) Remove IDENTICAL duplicate rows
%     TAB_tol_twin = unique(TAB_tol_twin, 'rows', 'stable');
% else
%     % first simulation
%     TAB_tol_twin = TAB;
% end
% 
% % 4) Sort (put here the columns that define the ordering)
% % Example: if the first column is an important parameter
% TAB_tol_twin = sortrows(TAB_tol_twin, TAB_tol_twin.Properties.VariableNames);
% 
% % 5) Always save to the same file
% save(fileName, 'TAB_tol_twin');
% % 
% % if isfile(fileName)
% %     delete(fileName)
% %     disp('TAB_tol_hoop.mat cancellato');
% % else
% %     disp('Il file non esiste');
% % end