function VIZfilenames = writeVizualizationFiles(simulationFolder, OP_file_or_struct, varargin)
% Loop through simulations, looking for *.postmb files, if present, VIZ files will be generated based on "best guess values", which may be overriden using varargin

% INPUTS:
%  - simulationFolder:  folder where linearization files are present.
%  - OP_file_or_struct: 
%          path to a csv file that contains information about the Operating points see function readOperatingPoints for more.
%       OR structure with (depending on simulation) fields: RotorSpeed, {optional: WindSpeed, GeneratorTorque, BladePitch, TowerTopDispFA, Filename}
%
% OPTIONAL INPUTS:
%  - varargin is a set of ('key', value) pairs used to override values found in VIZ file
%            Look at `opts` below for the different options, and default values.
%
% OUTPUTS:
%  - VIZfilenames: cell-array of VIZ filenames writen

% --- Optional arguments
OptsFields={'VTKLinModes','VTKModes','VTKLinScale','VTKLinTim','VTKLinTimes1','VTKLinPhase'};
opts=struct();
% Default values
opts.VTKLinModes = 15  ;% Number of modes to visualize (0 <= VTKLinModes <= NumModes)
opts.VTKModes    = '1,';% List of which VTKLinModes modes will be visualized (modes will be added sequentially from the last value entered)
opts.VTKLinScale = NaN; % Mode shape visualization scaling factor (exaggerates mode shapes). If NaN, the following is used: 10 for ElastoDyn; 0.1 for BeamDyn.
opts.VTKLinTim   = NaN; % Switch to make one animation for all LinTimes together (VTKLinTim=1) or separate animations for each LinTimes (VTKLinTim=2); If NaN, VTKLinTimes set to 2
opts.VTKLinTimes1= '' ; % If VTKLinTim=2, visualize modes at LinTimes(1) only? (if false, files will be generated at all LinTimes); If empty, set to true
opts.VTKLinPhase = 0  ; % If VTKLinTim=2, visualize modes at LinTimes(1) only? (if false, files will be generated at all LinTimes)
% Values input by users % NOTE: inputParser not available in Octave
if mod(length(varargin),2)~=0; error('Varargin should have an even number of values, for key/values pairs'); end
for iVar = 1:2:length(varargin)
    i = find( strcmp( OptsFields, varargin{iVar}) == 1);
    if ~isempty(i)
        opts.(OptsFields{i}) = varargin{iVar + 1};
    else
        warning('Optional key `%s` not supported by function %s',varargin{iVar},mfilename)
    end
end


% Getting fullpath to FASTfilenames
[FASTfilenames, OP] = getFullFilenamesOP(simulationFolder, OP_file_or_struct);


% Read VIZ template
%paramVIZ = FAST2Matlab(templateFilenameVIZ,2);


% Loop through simulations, looking for *.postmb files, if present FAST will be rerun to generate VTK
FileError={};
VIZfilenames ={};
for iOP = 1:length(FASTfilenames)
    fst = FASTfilenames{iOP};
    if isfield(OP,'RotorSpeed')
        RPM = OP.RotorSpeed(iOP);
    else
        RPM = 0;
    end

    % --- Read FAST file 
    FP = FAST2Matlab(fst,2); % Fast Parameters, specify 2 lines of header (FAST 8)
    CompElast  = GetFASTPar(FP,'CompElast') ;

    % --- Write VIZ file
    [fdir, base,~] = fileparts(fst);
    if isempty(fdir)
        fdir='.';
    end
    fullbase = strrep(strrep([fdir '/'  base], '//','/'),'\','/');
    full_chkpFile      = [fullbase '.ModeShapeVTK.chkp']   ; % Generated by FAST
    full_chkpRoot      = [fullbase '.ModeShapeVTK']        ; % No control over this root
    chkpRoot           = [base     '.ModeShapeVTK']        ; % No control over this root
    postMBC            = [base     '.ModeShapeVTK.postmbc']; 
    full_postMBC       = [fullbase '.ModeShapeVTK.postmbc']; % <<< NOTE: name set and generated by getCampbellData
    full_ModesVizName  = [fullbase '.ModeShapeVTK.viz']    ; % We will create that file

    if ~exist(full_postMBC, 'file')
        fprintf('Missing:    %s\n',full_postMBC);
    else
        % --- Patch template
        %paramVIZ_mod = SetFASTPar(paramVIZ    ,'CheckpointRoot',normpath(full_chkpRoot));
        %paramVIZ_mod = SetFASTPar(paramVIZ_mod,'MatlabFileName',normpath(full_postMBC ));
        % Override input file with user arguments
        %for iField = 1:length(OptsFields)
        %    field=OptsFields{iField};
        %    if ~isnan(opts.(field))
        %        paramVIZ_mod = SetFASTPar(paramVIZ_mod,field,opts.(field));
        %    end
        %end
        % TODO Check that VTKLinScale is about 10 foe CompElast=1 and around 0.1 for CompElast=2
        %paramVIZ_mod 
%             Matlab2FAST(paramVIZ_mod, templateFilenameVIZ, full_ModesVizName, 2); %contains 2 header lines
        % --- Write VIZ file directly
        try
            fid=fopen(full_ModesVizName, 'w');
        catch
            FileError{end+1} = full_ModesVizName;
            continue; % do loop
        end
        VIZfilenames{end+1}= normpath(full_ModesVizName); % ...
        fprintf(fid,'------- OpenFAST MODE-SHAPE INPUT FILE -------------------------------------------\n');
        fprintf(fid,'# Options for visualizing mode shapes\n');
        fprintf(fid,'---------------------- FILE NAMES ----------------------------------------------\n');
        fprintf(fid,'"%s"   CheckpointRoot - Rootname of the checkpoint file written when OpenFAST generated the linearization files (without the ".chkp" extension)\n',chkpRoot);
        fprintf(fid,'"%s"   MatlabFileName - Name of the mode-shape file (with eigenvectors) written by Matlab ("ModesVizName" sent to function fx_mbc3())\n'          ,postMBC);
        fprintf(fid,'---------------------- VISUALIZATION OPTIONS -----------------------------------\n');
        fprintf(fid,'%d        VTKLinModes   - Number of modes to visualize (0 <= VTKLinModes <= NumModes)\n',opts.VTKLinModes);
        fprintf(fid,'%s        VTKModes      - List of which VTKLinModes modes will be visualized (modes will be added sequentially from the last value entered)\n',opts.VTKModes);
        if isnan(opts.VTKLinScale)
            % Then user didn't specify it, we use some logic
            if CompElast==1 % ElastoDyn - VTKLinScale=10
                fprintf(fid,'10    VTKLinScale   - Mode shape visualization scaling factor (exaggerates mode shapes: try 10 for ElastoDyn; 0.1 for BeamDyn)\n');
            elseif CompElast==2 % BeamDyn - VTKLinScale=0.1
                fprintf(fid,'0.1   VTKLinScale   - Mode shape visualization scaling factor (exaggerates mode shapes: try 10 for ElastoDyn; 0.1 for BeamDyn)\n'); 
            end
        else
            if CompElast==1 && opts.VTKLinScale<1
                warning('VTKLinScale might be too small')
            end
            if CompElast==2 && opts.VTKLinScale>3
                warning('VTKLinScale might be too large')
            end
                fprintf(fid,'%f    VTKLinScale   - Mode shape visualization scaling factor (exaggerates mode shapes: try 10 for ElastoDyn; 0.1 for BeamDyn)\n',opts.VTKLinScale); 
        end
        if isnan(opts.VTKLinTim)
            % The user didn't specify this, we use some logic
            fprintf(fid,'2         VTKLinTim     - Switch to make one animation for all LinTimes together (VTKLinTim=1) or separate animations for each LinTimes (VTKLinTim=2)\n');
        else
            if (RPM<1e-3) % When RPM =0, VTKLinTim=1 would only produce one VTK
                fprintf(fid,'2         VTKLinTim     - Switch to make one animation for all LinTimes together (VTKLinTim=1) or separate animations for each LinTimes (VTKLinTim=2)\n');
            else
                fprintf(fid,'%d        VTKLinTim     - Switch to make one animation for all LinTimes together (VTKLinTim=1) or separate animations for each LinTimes (VTKLinTim=2)\n',opts.VTKLinTim);
            end
        end
        if length(opts.VTKLinTimes1)==0
            % The user didn't specify this
            fprintf(fid,'true      VTKLinTimes1  - If VTKLinTim=2, visualize modes at LinTimes(1) only? (if false, files will be generated at all LinTimes)\n');
        else
            fprintf(fid,'%s        VTKLinTimes1  - If VTKLinTim=2, visualize modes at LinTimes(1) only? (if false, files will be generated at all LinTimes)\n',opts.VTKLinTimes1);
        end
        fprintf(fid,'%f       VTKLinPhase   - Phase used when making one animation for all LinTimes together (used only when VTKLinTim=1)\n',opts.VTKLinPhase);
        fclose(fid);
        fprintf('Written:    %s\n',full_ModesVizName);
    end
end
nWritten = length(VIZfilenames)          ;
nError   = length(FileError)             ;
nMissing = length(FASTfilenames)-nWritten-nError;
fprintf('Viz files:  %d written, %d missing postmbc, %d failed\n',nWritten, nMissing, nError);
if length(FileError)>0
    disp('Error: the files below couldnt be written:')
    disp(FileError)
    error('Some VIZ files couldnt be written (see above)')
end

end

function s= normpath(s)
    if ispc()
        s = strrep(s,'/','\');
    else
        s = strrep(s,'\','/');
    end
end

