clc; clear; close all

% PARAMETRIC STUDY FOR SPRING-TYPE WHEEL (SPIRAL MODIFIED)

% GEOMETRIC PARAMETERS

g = 1.62;

%%% Values of the best configuration
R_ext_vec= 0.5;
B_vec= 0.4;
t_ext_vec= 0.0015 ;
t_int_vec= 0.002 ; 
t_rinf_vec= 0.002; 
N_vec= 8; 
M_vec= 7000; 

%%% Values for the Comparison analysis
% 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];

% t_ext_vec= [ 0.003 ]; % 0.008 0.012]; % outer ring thickness
% t_int_vec= [0.003]; % 0.008 0.012]; % inner ring thickness
% t_rinf_vec= [0.003]; % 0.004 0.006]; % reinforcement thickness 


% ARCHITECTURAL PARAMETERS (HOOP WHEEL)

R_int_vec= 0.25; % R_ext_vec/2;
R_media= R_int_vec+(R_ext_vec-R_int_vec)/2;
Circonf_vec= 2*pi*R_media;
N_elem_tot_vec= floor(Circonf_vec/(2*(R_ext_vec-R_int_vec)/2));

%N_elem_vec= floor(3:3:N_elem_tot_vec);
 N_elem_vec= 9;              % number of reinforcements

% 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= 0.025;

% 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);
                                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);

                                % SESSION FILE CREATION

                                file_ses= "architecture_tweel.ses";
                                file_sesw=fopen(file_ses,"w");
                                f_db_jou= "architecture_tweel.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_tweel.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',R_ext, ...
                                    '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, ...
                                    'mesh',mesh);
                                
                                
                                % PROPERTY STRUCTURE
                                
                                % Shell
                                param_Shell.name = "Shell_" + string(mat)';
                                param_Shell.th.name = ["t_ext", "t_int", "t_rinf"];
                                param_Shell.th.val= [t_ext,t_int,t_rinf];

                                % SESSION FILE PARAMETERS

                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, '$ GEOMETRY PARAMETERS\n');
                                fprintf(file_sesw, newline);
                                fprintf(file_sesw, 'REAL R_ext = %d\n', param_str.R_ext);
                                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);

                                % outer part

                                fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');

                                c_id=0;
                                
                                n_an= mod((0:360/param_str.N_elem:360)+270,360);  % point indices 
                                
                                for i=1:param_str.N_elem
                                    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_ext,n_an(i),n_an(i+1));
                                    c_id=c_id+1;
                                end
                                
                                p_id=c_id;
                                list_p_non_extr=1: p_id;
                                
                                
                                fprintf(file_sesw, 'STRING sgm_sweep_surface_e_created_ids[VIRTUAL]\n');
                                
                                fprintf(file_sesw, 'sgm_const_surface_extrude( "1", "<0 0 %d>", 1., 0., "[0 0 0]", "Coord 0", "Curve 1:%i", sgm_sweep_surface_e_created_ids )\n',param_str.B,c_id);
                                s_id=c_id;
                                p_id=p_id+s_id;
                                list_p_ext=list_p_non_extr(end)+1:p_id;
                                
                                for i=1:param_str.N_elem
                                    fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, %d, %d, "Coord 0.3", "[0 0 %d]", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_ext,n_an(i),n_an(i+1),param_str.B);
                                    c_id=c_id+1;
                                end
                                
                                list_p=1:p_id;
                                list_c=1:c_id;
                                list_s=1:s_id;

                                % hub
                                
                                for i=1:param_str.N_elem
                                    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+1));
                                    c_id=c_id+1;
                                    p_id=p_id+1;
                                end
                                list_p_int_non_extr= list_p(end)+1:p_id; 
                                
                                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,list_c(end)+1,c_id);
                                s_id=s_id+param_str.N_elem;
                                
                                for i=1:param_str.N_elem
                                    fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, %d, %d, "Coord 0.3", "[0 0 %d]", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1,param_str.R_int,n_an(i),n_an(i+1),param_str.B);
                                    c_id=c_id+1;
                                    p_id=p_id+1;
                                end
                                list_p_int_extr= list_p_int_non_extr(end)+1:p_id;
                                
                                list_p_int= list_p(end)+1:p_id;
                                list_c_int=list_c(end)+1:c_id;
                                list_s_int=list_s(end)+1:s_id;


                                % create reinforcements 
                                % create points for reinforcements
                                
                                fprintf(file_sesw, newline);

                                for i=1:param_str.N_elem

                                    if i==1

                                        fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, 40., 400., "Coord 0.3", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1, param_str.R_ext/10,list_p_int_non_extr(i));
                                        c_id=c_id+1;  % creates one extra curve for each element
                                        p_id=p_id+1;  % creates one extra point for each element that I do not need!
                                            
                                        fprintf(file_sesw, 'STRING asm_create_grid_int_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'asm_const_grid_pierce( "%i", "Curve %i", "Surface %i %i", asm_create_grid_int_created_ids )\n',p_id+1,c_id,list_s_int(1),list_s_int(end));
                                        p_id=p_id+2;

                                    elseif i>1 

                                        fprintf(file_sesw, 'STRING sgm_create_curve_2d_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'sgm_const_curve_2d_arcangles_v1( "%i", %d, 40., 400., "Coord 0.3", "Point %i", FALSE, sgm_create_curve_2d_created_ids )\n',c_id+1, param_str.R_ext/10,list_p_int_non_extr(i));
                                        c_id=c_id+1;  % creates one extra curve for each element
                                        p_id=p_id+1;  % creates one extra point for each element that I do not need!

                                        fprintf(file_sesw, 'STRING asm_create_grid_int_created_ids[VIRTUAL]\n');
                                        fprintf(file_sesw, 'asm_const_grid_pierce( "%i", "Curve %i", "Surface %i %i", asm_create_grid_int_created_ids )\n',p_id+1,c_id,list_s_int(i-1),list_s_int(i));
                                        p_id=p_id+2;

                                    end

                                end

                                list_p_rinf = list_p_int(end)+1:p_id;  % total point list (also including those useless n_elem points)
                                list_c_rinf_inutili= list_c_int(end)+1:c_id; % list of curves needed only to create points ... then I will delete them

                                list_p_rinf_inutili= list_p_rinf(1):3:list_p_rinf(end); % list of points that I will delete
                                list_p_rinf_utili= setdiff(list_p_rinf,list_p_rinf_inutili); % list of points that I need 

                                fprintf(file_sesw, 'STRING asm_delete_any_deleted_ids[VIRTUAL]\n');   % delete useless curves and points
                                fprintf(file_sesw, 'asm_delete_curve( "Curve %i:%i", asm_delete_any_deleted_ids )\n',list_c_rinf_inutili(1),list_c_rinf_inutili(end));

                                p_id= p_id-param_str.N_elem; % number of actually existing points (i.e. without the deleted ones)
                                c_id= c_id - length(list_c_rinf_inutili); % number of actually existing curves 

                                % create reinforcement curves 

                                for i=1:param_str.N_elem

                                    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_p_non_extr(i),list_p_rinf_utili(i));
                                    c_id=c_id+1;
                                    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_p_non_extr(i),list_p_rinf_utili(i+1));
                                    c_id=c_id+1;
                                    list_p_rinf_utili= setdiff(list_p_rinf_utili,list_p_rinf_utili(1));

                                end

                                list_p_rinf_utili= setdiff(list_p_rinf,list_p_rinf_inutili);
                                list_c_rinf = list_c_int(end)+1:c_id;
                                
                                % create surfaces for reinforcements
                                
                                for i=1:param_str.N_elem*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", sgm_sweep_surface_e_created_ids )\n',s_id+1,param_str.B,list_c_rinf(i));
                                    s_id=s_id+1;

                                end

                                list_s_rinf=list_s_int(end)+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+1,param_str.B/2);
                                p_id=p_id+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_ext=Al, t_int=Ac, t_rinf=Pol

                                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_int, list_s_rinf};   % surface lists

                                for k=1:n_th

                                    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);

                                fprintf(file_sesw, 'loadsbcs_create2( "Incastro", "Displacement", "Nodal", "", "Static", ["Point %i %i"], "Geometry", "Coord 0", "1.", ["< 0 0 0 >", "<  0 0 0   >", "<     >", "<     >"], ["", "", "", ""] )\n', list_p(1),list_p(1+N_elem));

                                fprintf(file_sesw, 'loadsbcs_create2( "Force_g", "Force", "Nodal", "", "Static", ["Point %i"], "Geometry", "Coord 0", "1.", ["<    0, %d, 0 >", "<     >", "<     >", "<     >"], ["", "", "", ""] )\n', p_id,-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);

                                % 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_int(1),list_s_int(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_tweel", "ANALYZE NO JOBFILE" )\n');
                                fprintf(file_sesw, 'msc_delete_old_files( "architecture_tweel", ".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_tweel\\architecture_tweel.db" )\n');
                                fprintf(file_sesw, 'jobfile.writec( "JOBNAME", "architecture_tweel" )\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_tweel", 1, ["Default"] )\n');
                                fprintf(file_sesw, 'analysis_submit_2( "MSC.Nastran", "architecture_tweel" )\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_tweel.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_tweel.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; 
                                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(3);
                                MATRIX(count, 12) = U_Y;
                                MATRIX(count, 13) = S_Mises_max;
                                MATRIX(count, 14) = mass;
                                MATRIX(count, 15) = K;
                                MATRIX(count, 16) = DI;
                                MATRIX(count, 17) = mesh;

                                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', 'Mesh'};

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


%%% TAB for the sensitivity analysis
% folderName = 'C:\Users\Sofia\OneDrive - Politecnico di Torino\TESI MAGISTRALE\Matlab tesi\Cartella risultati analisi';
% TAB_all_tweel_new= TAB;
% fileName = fullfile(folderName, 'TAB_all_tweel_new.mat');
% 
% save(fileName, 'TAB_all_tweel_new')


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

% if isfile(fileName)
%     delete(fileName)
%     disp('TAB_tol_hoop.mat cancellato');
% else
%     disp('Il file non esiste');
% end

% writetable(TAB_tweel_final,'TAB_tweel_final.xlsx');


%%% Case 1: TAB for the convergence analysis (best case)

% folderName = 'C:\Users\Sofia\OneDrive - Politecnico di Torino\TESI MAGISTRALE\Matlab tesi\Cartella risultati analisi\tweel_risultati';
% 
% fileName = fullfile(folderName, 'TAB_tweel_convergence.mat');
% 
% % 1) If it exists, load the old table
% if isfile(fileName)
%     S = load(fileName, 'TAB_tweel_convergence');
%     TAB_tweel_convergence = S.TAB_tweel_convergence;
% 
%     % 2) Append the new simulation
%     TAB_tweel_convergence = [TAB_tweel_convergence; TAB];
% 
%     % 3) Remove IDENTICAL duplicate rows
%     TAB_tweel_convergence = unique(TAB_tweel_convergence, 'rows', 'stable');
% else
%     % first simulation
%     TAB_tweel_convergence = TAB;
% end

% 4) Sort (put here the columns that define the ordering)
% Example: if the first column is an important parameter
% TAB_tweel_convergence = sortrows(TAB_tweel_convergence, TAB_tweel_convergence.Properties.VariableNames);

% % 5) Always save to the same file
% save(fileName, 'TAB_tweel_convergence');

%%% Case 2: Limit case

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

%%% Case 3 (lighter)

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