text:
% Write Values to the Windows Registry:
%  writereg(rootkey,subkey,valname,value,...)
%
%  Example:
%  writereg('HKEY_LOCAL_MACHINE','SOFTWARE\My Software, Inc.\AppName',...
%           'TextSetting','c:\0.txt',...
%           'NumSetting',1);
%
%  The following registry value types are used:
%    "REG_SZ": when value is a string.
%    "REG_DWORD": when value is an integer.
%  No other value types are supported.
%
%  See also: winqueryreg
%
%  Alex Nelson, August 29, 2007
%
%  This function is available without warranty. It is recommended that
%  you comment out the last line and verify the correctness of the .reg
%  file before use, and that you backup your windows registry.
%
%  This function was based on code available at:
%  www.compactech.com/forum/index.php?showtopic=241&mode=threaded&pid=590
%
function writereg(rootkey,subkey,varargin)
  if mod(length(varargin),2)
    error('writereg:Value names and values must be in pairs');
  end
  tmpregfile = [getenv('TEMP') '\write_reg.reg'];
  fp = fopen(tmpregfile,'wt');
  fprintf(fp,'Windows Registry Editor Version 5.00\n\n');
  fprintf(fp,'[%s\\%s]\n',rootkey,subkey);

  ii = 1;
  while ii<length(varargin)
    valname =varargin{ii};
    value = varargin{ii+1};
    if ischar(value)                       % string value
      value = strrep(value,'\','\\');     %  escape backslashes
      vstr=sprintf('"%s"="%s"',valname,value);
    elseif isint(value) && value<2^32     % 32-bit integer values
      vstr=sprintf('"%s"=dword:%.8X',valname,value); % write as hexidecimal
    else
      error('Only values of type string or integer are supported');
    end
    ii=ii+2;
    fprintf(fp,'%s\n',vstr);
  end
  fprintf(fp,'\n');
  fclose(fp);

  doscmd = sprintf('C:\\windows\\regedit.exe /s %s', tmpregfile)
  dos(doscmd);

% ISINT(A) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Returns 1 for entries in matrix which are integers, and 0 for entries
%  which are non-integer.
% Alex Nelson 12/18/03
function a = isint(A)
  a = A==fix(A);

author: Sable

repo: mcbench-benchmarks


text:
function modelnet_binlist2hdf5(DestinationFolder,Pose_plan)

Pose_file=['../poseplans/' Pose_plan '.txt'];
if ~exist(strcat(DestinationFolder,Pose_plan,'/'), 'dir')
    mkdir(strcat(DestinationFolder,Pose_plan,'/'));
end
copyfile(Pose_file,strcat(DestinationFolder,Pose_plan,'/'));

%% ---- Create the lists
rotation_start_index = 1; %this is the original Shapenet convention
collate_singlerots = 0; %it has no effect for modelnet10 , with the current poseplan
alllist_file = sprintf('%s/%s/all.txt',DestinationFolder,Pose_plan);
system(sprintf('find %s -iname ''*.%s'' | sort -V > %s', DestinationFolder,'bin',alllist_file));
modelnet_add_labels_to_list(alllist_file,alllist_file,rotation_start_index,collate_singlerots,Pose_file);

lists_dir = sprintf('%s/%s/lists',DestinationFolder,Pose_plan);
prepare_modelnet_lists(alllist_file,lists_dir,rotation_start_index);

%% --- Create some special HDF5 datasets
binlist_to_hdf5([lists_dir '/shuffled/train_allrot.txt'],[DestinationFolder '/' Pose_plan '/hdf5/train_allrot_shuffled/train.hdf5'],[DestinationFolder '/' Pose_plan '/hdf5/train_allrot_shuffled/train.hdf5.txt'],1000);
binlist_to_hdf5([lists_dir '/test_singlerandrot.txt'],[DestinationFolder '/' Pose_plan '/hdf5/test_singlerandrot/test.hdf5'],[DestinationFolder '/' Pose_plan '/hdf5/test_singlerandrot/test.hdf5.txt'],1000);
binlist_to_hdf5([lists_dir '/test_allrot.txt'],[DestinationFolder '/' Pose_plan '/hdf5/test_allrot/test.hdf5'],[DestinationFolder '/' Pose_plan '/hdf5/test_allrot/test.hdf5.txt'],1000);

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
%% (Internal) Write annotation files for biomedical signals in MIT Format. (MEX file)
%
% Description
% Input ann_struct is a struct with 6 fields: time, anntyp, subtyp, chan, num
% and aux, as the one returned by readannot.
% Each of this fields is a vector with as many elements as the number of events
% in the anotator.  The struct can be incomplete, but time and anntyp fields are
% strictly necessary.
%
%     writeannot(ann_filename,ann_struct)
%
% Arguments:
%
%   + ECG_header: original header struct.
%
%   + ECG_idx: signal indexes to select.
%
% Output:
%
%   + ECG_header: trimed header struct.
%
% See also read_ecg, writeheader, ECGwrapper
%
% Author: Juan Pablo Martinez Cortes
% adapted by Mariano Llamedo Soria
% Version: 0.1 beta
% Birthdate: 21/7/2010
% Last update: 20/02/2013
% Copyright 2008-2015
%

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [propF, propTS, logQ_RGS] = sampleFeatures_RGS( F, TS, theta, data_struct, hyperparams, model, anchorObjIDs, propFeatIDs, TargetPsi )
% Sample binary feature asgn matrix F  under RESTRICTED settings
%   which means we
%       (1) only update features for sequences in list "objIDs"
%       (2) only consider asgns where at least one of "featIDs" is positive
% Executed under two circumstances:
%   (I)  actually sample from posterior on F over [0 1], [1 0], [1 1]
%   (II) do not sample at all, but compute prob.
%            for moving from current F to TargetPsi.F with fixed TS,theta

propF = F;
logQ_RGS = 0;
obsModel = model.obsModel;

activeObjIDs = find( sum(F(:, propFeatIDs), 2) > 0  )';

% Expand TS to make sure all objects possess both feature ids
propTS = expandTStoIncludeFeatIDs( F, TS, hyperparams, activeObjIDs, propFeatIDs );

% Initialize some key sufficient stats
featCounts = sum( F( activeObjIDs, propFeatIDs ), 1 );
nObj = size(F,1);

didWarn=0;
for aa = activeObjIDs

     % Make sure to update F suff. stats
     featCounts = featCounts - propF(aa, propFeatIDs);
     nObj = nObj - 1;

    % Calc Soft Evidence for the current sequence, including both new feats
    curFeatIDs = find( propF(aa,:));
    allPropFeatIDs = union( curFeatIDs, propFeatIDs );
    LL = calcLogCondPrObsGivenTheta( data_struct(aa), theta, obsModel, allPropFeatIDs, size(propF,2) );

    % Propose new feature assignments for the cur. seq.
    myFPr = featCounts ./ ( nObj + hyperparams.c );
    [newFaa, logQ_Faa] = sampleFfromRestrictedSet( aa, anchorObjIDs, LL, propFeatIDs, propTS.obj( aa ), myFPr, TargetPsi );
    propF(aa, propFeatIDs) = newFaa;

    % Keep track of probability of making this specific proposal
    logQ_RGS = logQ_RGS + logQ_Faa;
    if ~isinf( logQ_RGS ) && ~didWarn
        %fprintf( 'WARNING: Bad Calculation of Proposal Prob. for F\n' );
        didWarn=1;
    end
    % In case we chose not to include one of the proposed new split feats,
    %   we need to remove that feature's transition weights from propTS
    %   and also remove that feature's entries in soft ev. matrix LL
    availFeatIDs = propTS.obj(aa).availFeatIDs;
    kiloc = find( availFeatIDs == propFeatIDs(1) );
    kjloc = find( availFeatIDs == propFeatIDs(2) );
    switch ( newFaa * [10;1] )
        case 01
            keepIDs = [1:kiloc-1 kiloc+1:length(availFeatIDs)];
        case 10
            keepIDs = [1:kjloc-1 kjloc+1:length(availFeatIDs)];
        case 11
            keepIDs = 1:length(availFeatIDs);
    end
    propTS.obj(aa).availFeatIDs = propTS.obj(aa).availFeatIDs( keepIDs );
    propTS.obj(aa).pi_z = propTS.obj(aa).pi_z( keepIDs, keepIDs );
    propTS.obj(aa).pi_init = propTS.obj(aa).pi_init( keepIDs );

    % Update suff. stats for F sampling
    featCounts = featCounts + propF(aa, propFeatIDs);
    nObj = nObj + 1;

end

author: michaelchughes

repo: NPBayesHMM


text:
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
figure;
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
clear mycolor
clear all_clusters
clear all_clusterHandles

%------------------------------------------------------------------------------
% Get head, if not default
%------------------------------------------------------------------------------
cflag = 0;
while cflag == 0
    a = input('Use default head single_subj_T1...?  (type default or file name, no .img extension)','s');
    switch a
    case 'default', fileName = 'single_subj_T1';
    otherwise fileName = a;
    end
    if exist([fileName '.img']) == 2,
	cflag = 1;
    else
	disp(['Can''t locate file: ' fileName '.img'])
    end
end

%------------------------------------------------------------------------------
% Add surface?
%------------------------------------------------------------------------------
cflag = 0;
while cflag == 0
    a = input('Add brain surface (y\n\filename,no .img extension)...?','s');
    switch a
    case 'y', fileName2 = ['surf_' fileName];, dosurf = 1;
    case 'n', fileName2 = [];, cflag = 1;, dosurf = 0;
    otherwise fileName2 = a;, dosurf = 1;
    end
    if dosurf & exist([fileName2 '.mat']) == 2,
	    cflag = 1;
    elseif dosurf
	    disp(['Can''t locate file: ' fileName '.mat'])
    end
end

numSets = input('Enter number of cluster sets: ');

for i = 1:numSets

    disp(['Cluster Set ' num2str(i)])
    disp('_______________________________________________')

    %------------------------------------------------------------------------------
    % Get clusters, if necessary
    %------------------------------------------------------------------------------
    cflag = 0;
    while cflag == 0
       a = input('Get clusters from...?  (workspace,file,figure)','s');
       switch a
       case 'workspace',
            cflag = 1;
            varname = input('Enter name of workspace variable: ','s');
            eval(['clusters = ' varname ';'])
            if isempty(clusters),disp('No clusters or not a valid variable name.'),cflag=0;,end
       case 'file', clusters = imageCluster('getclusters');, cflag = 1;
       case 'figure', clusters = imageCluster('getfigclusters');, cflag = 1;
       otherwise disp(['Not a valid choice...type entry string with no quotes.'])
       end
    end

    mycolor{i} = input('Enter color string in single quotes or color vector for clusters: ');
    all_clusters{i} = clusters;
end

%------------------------------------------------------------------------------
% Make 3 views; loop through cluster sets
%------------------------------------------------------------------------------

figure;
myView = {[0 90] [270 0] [180 0]}

for myAxis = 1:3
    subplot(1,3,myAxis); hold on

    for i = 1:numSets
        cH = [];
        for j = 1:length(all_clusters{i})
            cH(j) = imageCluster('cluster',all_clusters{i}(j),'color',mycolor{i},'alpha',.6); drawnow
            set(cH(j),'Tag',['set' num2str(i) 'cluster' num2str(j)])
        end

        clusterHandles{i} = cH;

    end

    all_clusterHandles{myAxis} = clusterHandles;

    %------------------------------------------------------------------------------
    % * lighting stuff
    %------------------------------------------------------------------------------
    view(myView{myAxis}(1),myView{myAxis}(2))
    lighting gouraud;
    [az,el] = view;
    myLight = lightangle(az,el);
    set(myLight,'Tag','myLight')
    rotate3d on

    % set callback to light follow the camera
    %------------------------------------------------------------------------------
    if exist('lightFollowView') == 2
         set(gcf, 'WindowButtonUpFcn', 'lightFollowView');
    else
        warning('Cannot find lightFollowView.m to set light position.')
    end
    axis vis3d
    axis image
    drawnow

    %------------------------------------------------------------------------------
    % * head surface
    %------------------------------------------------------------------------------
    if dosurf
        p = [];
        eval(['load ' fileName2])
        try
            p(end+1) = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5],'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',.2,'SpecularExponent',200);
        catch
            p(end+1) = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5],'EdgeColor','none','SpecularStrength',.2,'SpecularExponent',200);
        end
        axis vis3d
        axis image
        drawnow
    end

end     % loop through axes

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
report_this_filefun(mfilename('fullpath'));

X = reshape(c(:,1),281,221);
Y = reshape(c(:,2),281,221);
Z = reshape(c(:,3),281,221);

figure
pcolor(X,Y,(Z));shading interp

author: CelsoReyes


text:
%%*****************************************************************************
%% HSDsqlp: solve an semidefinite-quadratic-linear program
%%   by infeasible path-following method on the homogeneous self-dual model.
%%
%%  [obj,X,y,Z,info,runhist] =
%%      HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0,kap0,tau0,theta0);
%%
%%  Input:
%%  blk    : a cell array describing the block diagonal structure of SQL data.
%%   At    : a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%  b,C    : data for the SQL instance.
%%  OPTIONS: a structure that specifies parameters required in HSDsqlp.m,
%%           (if it is not given, the default in sqlparameters.m is used).
%%
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%  (kap0,tau0,theta0): initial parameters (if not given, the default is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%----------------------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  ...
%%      (all have default values set in sqlparameters.m).
%%
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

function [obj,X,y,Z,info,runhist] = ...
    HSDsqlpmain(blk,At,C,b,par,X0,y0,Z0,kap0,tau0,theta0)

%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%
global spdensity  smallblkdim  solve_ok use_LU
global schurfun   schurfun_par
% matlabversion = par.matlabversion;
isoctave = exist( 'OCTAVE_VERSION', 'builtin' );
if isoctave,
  w1 = warning('off','Octave:nearly-singular-matrix');
else
  w1 = warning('off','MATLAB:nearlySingularMatrix');
  w2 = warning('off','MATLAB:singularMatrix');
end
vers          = par.vers;
predcorr      = par.predcorr;
gam           = par.gam;
expon         = par.expon;
gaptol        = par.gaptol;
inftol        = par.inftol;
steptol       = par.steptol;
maxit         = par.maxit;
printlevel    = par.printlevel;
stoplevel     = par.stoplevel;
% scale_data    = par.scale_data;
spdensity     = par.spdensity;
rmdepconstr   = par.rmdepconstr;
smallblkdim   = par.smallblkdim;
schurfun      = par.schurfun;
schurfun_par  = par.schurfun_par;
% ublksize      = par.ublksize;
%%
tstart = clock;
X = X0; y = y0; Z = Z0;
for p = 1:size(blk,1)
    if strcmp(blk{p,1},'u'); Z{p} = zeros(blk{p,2},1); end
end
%%
%%-----------------------------------------
%% convert unrestricted blk to linear blk.
%%-----------------------------------------
%%
ublkidx = zeros(size(blk,1),1);
Cpert = zeros(size(blk,1),1);
Cnew = C;
perturb_C = 1;
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    tmp = max(1,norm(C{p},'fro'))/sqrt(n);
    if strcmp(pblk{1},'s')
        if (perturb_C); Cpert(p) = 1e-3*tmp; end
        Cnew{p} = C{p} + Cpert(p)*speye(n);
    elseif strcmp(pblk{1},'q')
        if (perturb_C); Cpert(p) = 0*tmp; end; %% old: 1e-3
        s = 1+[0, cumsum(pblk{2})];
        tmp2 = zeros(n,1); len = length(pblk{2});
        tmp2(s(1:len)) = ones(len,1);
        Cnew{p} = C{p} + Cpert(p)*tmp2;
    elseif strcmp(pblk{1},'l')
        if (perturb_C); Cpert(p) = 1e-4*tmp; end; %% old: 1e-3
        Cnew{p} = C{p} + Cpert(p)*ones(n,1);
    elseif strcmp(pblk{1},'u')
        msg = sprintf('convert ublk to linear blk');
        if (printlevel); fprintf('\n *** %s',msg); end
        ublkidx(p) = 1;
        n = 2*pblk{2};
        blk{p,1} = 'l'; blk{p,2} = n;
        if (perturb_C); Cpert(p) = 1e-2*tmp; end
        C{p}  = [C{p}; -C{p}];
        At{p} = [At{p}; -At{p}];
        Cnew{p} = C{p} + Cpert(p)*ones(n,1);
        X{p}  = 1+randmat(n,1,0,'u');
        Z{p}  = 1+randmat(n,1,0,'u');
    end
end
%%
%%-----------------------------------------
%% check if the matrices Ak are
%% linearly independent.
%%-----------------------------------------
%%
m0 = length(b);
[At,b,y,indeprows,depconstr,feasible,AAt] = ...
    checkdepconstr(blk,At,b,y,rmdepconstr);
if (~feasible)
    msg = 'SQLP is not feasible';
    if (printlevel); fprintf('\n %s',msg); end
    return;
end
par.depconstr = depconstr;
%%
normb = 1+max(abs(b));
normC = zeros(length(C),1);
for p = 1:length(C); normC(p) = max(max(abs(C{p}))); end
normC = 1+max(normC);
nn = ops(C,'getM');
m = length(b);
if (nargin <= 8) || (isempty(kap0) || isempty(tau0) || isempty(theta0))
    if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e6)
        kap0 = 10*blktrace(blk,X,Z);
    else
        kap0 = blktrace(blk,X,Z);
    end
    tau0 = 1; theta0 = 1;
end
kap = kap0; tau = tau0; theta = theta0;
%%
normX0 = ops(X0,'norm')/tau; normZ0 = ops(Z0,'norm')/tau;
bbar  = (tau*b-AXfun(blk,At,[],X))/theta;
ZpATy = ops(Z,'+',Atyfun(blk,At,[],[],y));
Cbar = ops(ops(ops(tau,'*',C),'-',ZpATy),'/',theta);
gbar = (blktrace(blk,C,X)-b'*y+kap)/theta;
abar = (blktrace(blk,X,Z)+tau*kap)/theta;
for p = 1:size(blk,1);
    pblk = blk(p,:);
    if strcmp(pblk{1},'s')
        At{p} = [At{p}, -svec(pblk,Cnew{p},1), svec(pblk,Cbar{p},1)];
    else
        At{p} = [At{p}, -Cnew{p}, Cbar{p}];
    end
end
Bmat = [sparse(m,m), -b, bbar; b', 0, gbar; -bbar', -gbar, 0];
em1 = zeros(m+2,1); em1(m+1) = 1;
em2 = zeros(m+2,1); em2(m+2) = 1;
par.Umat = [[b;0;0], [bbar;gbar;0], em1, em2];
par.m = m;
par.diagAAt = [full(diag(AAt)); 1; 1];
%%
%%-----------------------------------------
%% find the combined list of non-zero
%% elements of Aj, j = 1:k, for each k.
%%-----------------------------------------
%%
par.numcolAt = length(b)+2;
[At,C,Cnew,X,Z,par.permA,par.invpermA,par.permZ] = ...
    HSDsortA(blk,At,C,Cnew,[b;0;0],X,Z); %#ok
[par.isspA,par.nzlistA,par.nzlistAsum,par.isspAy,par.nzlistAy] = ...
    nzlist(blk,At,par);
%%
%%-----------------------------------------
%% initialization
%%-----------------------------------------
%%
y2 = [y; tau; theta];
AX = AXfun(blk,At,par.permA,X);
rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
% Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
trXZ = blktrace(blk,X,Z);
mu  = (trXZ+kap*tau)/(nn+1);
obj = [blktrace(blk,C,X), b'*y]/tau;
gap = trXZ/tau^2;
relgap = gap/(1+mean(abs(obj)));
ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,[y;0;0]));
ZpATynorm = ops(ZpATy,'norm');
prim_infeas = norm(b - AX(1:m)/tau)/normb;
dual_infeas = ops(ops(C,'-',ops(ZpATy,'/',tau)),'norm')/normC;
infeas = max(prim_infeas,dual_infeas);
%%
termcode = 0;
pstep = 1; dstep = 1; pred_convg_rate = 1; corr_convg_rate = 1;
besttol = max( relgap, infeas );
% homRd = inf; homrp = inf; dy = zeros(length(b),1);
msg = []; msg2 = [];
runhist.pobj    = obj(1);
runhist.dobj    = obj(2);
runhist.gap     = gap;
runhist.relgap  = relgap;
runhist.pinfeas = prim_infeas;
runhist.dinfeas = dual_infeas;
runhist.infeas  = infeas;
runhist.cputime = etime(clock,tstart);
runhist.step    = 0;
runhist.kappa   = kap;
runhist.tau     = tau;
runhist.theta   = theta;
runhist.useLU   = 0;
ttime.preproc   = runhist.cputime;
ttime.pred = 0; ttime.pred_pstep = 0; ttime.pred_dstep = 0;
ttime.corr = 0; ttime.corr_pstep = 0; ttime.corr_dstep = 0;
ttime.pchol = 0; ttime.dchol = 0; ttime.misc = 0;
%%
%%-----------------------------------------
%% display parameters, and initial info
%%-----------------------------------------
%%
if (printlevel >= 2)
    fprintf('\n********************************************');
    fprintf('************************************************\n');
    fprintf('   SDPT3: homogeneous self-dual path-following algorithms');
    fprintf('\n********************************************');
    fprintf('************************************************\n');
    [hh,mm,ss] = mytime(ttime.preproc);
    if (printlevel>=3)
        fprintf(' version  predcorr  gam  expon\n');
        if (vers == 1); fprintf('   HKM '); elseif (vers == 2); fprintf('    NT '); end
        fprintf('     %1.0f      %4.3f   %1.0f\n',predcorr,gam,expon);
        fprintf('it pstep dstep pinfeas dinfeas  gap')
        fprintf('     mean(obj)    cputime    kap   tau    theta\n');
        fprintf('------------------------------------------------');
        fprintf('--------------------------------------------\n');
        fprintf('%2.0f|%4.3f|%4.3f|%2.1e|%2.1e|',0,0,0,prim_infeas,dual_infeas);
        fprintf('%2.1e|%- 7.6e| %s:%s:%s|',gap,mean(obj),hh,mm,ss);
        fprintf('%2.1e|%2.1e|%2.1e|',kap,tau,theta);
    end
end
%%
%%---------------------------------------------------------------
%% start main loop
%%---------------------------------------------------------------
%%
EE = ops(blk,'identity');
normE = ops(EE,'norm');
Zpertold = 1;
[Xchol,indef(1)] = blkcholfun(blk,X);
[Zchol,indef(2)] = blkcholfun(blk,Z);
if any(indef)
    msg = 'stop: X, Z are not both positive definite';
    if (printlevel); fprintf('\n  %s\n',msg); end
    info.termcode = -3;
    info.msg1 = msg;
    return;
end
%%
param.termcode    = termcode;
param.iter        = 0;
param.normX0      = normX0;
param.normZ0      = normZ0;
param.m0          = m0;
param.indeprows   = indeprows;
param.prim_infeas_bad = 0;
param.dual_infeas_bad = 0;
param.prim_infeas_min = prim_infeas;
param.dual_infeas_min = dual_infeas;
param.gaptol      = gaptol;
param.inftol      = inftol;
param.maxit       = maxit;
param.printlevel  = printlevel;
param.stoplevel   = stoplevel;
breakyes = 0;
dy = zeros(length(b),1); dtau = 0; dtheta = 0;
Xbest = X; ybest = y; Zbest = Z;
% kapbest = kap; taubest = tau; thetabest = theta;
%%
for iter = 1:maxit;
    update_iter = 0; pred_slow = 0; corr_slow = 0; % step_short = 0;
    tstart  = clock;
    timeold = tstart;
    par.kap   = kap;
    par.tau   = tau;
    par.theta = theta;
    par.mu   = mu;
    par.iter = iter;
    par.y    = y;
    par.dy2  = [dy; dtau; dtheta];
    par.rp   = rp;
    par.ZpATynorm = ZpATynorm;
    %%
    %%--------------------------------------------------
    %% perturb C
    %%--------------------------------------------------
    %%
    if (perturb_C)
        [At,Cpert] = HSDsqlpCpert(blk,At,par,C,X,Cpert,runhist);
        maxCpert(iter) = max(Cpert); %#ok
        %%fprintf(' %2.1e',max(Cpert));
        if (iter > 10 && norm(diff(maxCpert([iter-3,iter]))) < 1e-13)
            Cpert = 0.5*Cpert;
            maxCpert(iter) = max(Cpert); %#ok
        end
        AX = AXfun(blk,At,par.permA,X);
        rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
        Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
    end
    %%---------------------------------------------------------------
    %% predictor step.
    %%---------------------------------------------------------------
    %%
    if (predcorr)
        sigma = 0;
    else
        sigma = 1-0.9*min(pstep,dstep);
        if (iter == 1); sigma = 0.5; end;
    end
    sigmu = sigma*mu;

    invXchol = cell(size(blk,1),1);
    invZchol = ops(Zchol,'inv');
    if (vers == 1);
        [par,dX,dy,dZ,coeff,L,hRd] = ...
            HSDHKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol);
    elseif (vers == 2);
        [par,dX,dy,dZ,coeff,L,hRd] = ...
            HSDNTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);
    end
    if (solve_ok <= 0)
        msg = 'stop: difficulty in computing predictor directions';
        if (printlevel); fprintf('\n  %s',msg); end
        runhist.cputime(iter+1) = etime(clock,tstart);
        termcode = -4;
        break;
    end
    timenew = clock;
    ttime.pred = ttime.pred + etime(timenew,timeold); timeold=timenew;
    %%
    %%-----------------------------------------
    %% step-lengths for predictor step
    %%-----------------------------------------
    %%
    if (gam == 0)
        gamused = 0.9 + 0.09*min(pstep,dstep);
    else
        gamused = gam;
    end
    kapstep = max( (par.dkap<0)*(-kap/(par.dkap-eps)), (par.dkap>=0)*1e6 );
    taustep = max( (par.dtau<0)*(-tau/(par.dtau-eps)), (par.dtau>=0)*1e6 );
    [Xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);
    timenew = clock;
    ttime.pred_pstep = ttime.pred_pstep + etime(timenew,timeold); timeold=timenew;
    Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
    pstep = min(1,gamused*min([Xstep,Zstep,kapstep,taustep]));
    dstep = pstep;
    kappred = kap + pstep*par.dkap;
    taupred = tau + pstep*par.dtau;
    trXZpred = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ) ...
        + pstep*dstep*blktrace(blk,dX,dZ);
    mupred  = (trXZpred + kappred*taupred)/(nn+1);
    mupredhist(iter) = mupred; %#ok
    timenew = clock;
    ttime.pred_dstep = ttime.pred_dstep + etime(timenew,timeold); timeold=timenew;
    %%
    %%-----------------------------------------
    %%  stopping criteria for predictor step.
    %%-----------------------------------------
    %%
    if (min(pstep,dstep) < steptol) && (stoplevel)
        msg = 'stop: steps in predictor too short';
        if (printlevel)
            fprintf('\n  %s',msg);
            fprintf(': pstep = %3.2e,  dstep = %3.2e',pstep,dstep);
        end
        runhist.cputime(iter+1) = etime(clock,tstart);
        termcode = -2;
        breakyes = 1;
    end
    if (iter >= 2)
        idx = max(2,iter-2) : iter;
        pred_slow = all(mupredhist(idx)./mupredhist(idx-1) > 0.4);
        idx = max(2,iter-5) : iter;
        pred_convg_rate = mean(mupredhist(idx)./mupredhist(idx-1));
        pred_slow = pred_slow + (mupred/mu > 5*pred_convg_rate);
    end
    if (~predcorr)
        if (max(mu,infeas) < 1e-6) && (pred_slow) && (stoplevel)
            msg = 'stop: lack of progress in predictor';
            if (printlevel)
                fprintf('\n  %s',msg);
                fprintf(': mupred/mu = %3.2f, pred_convg_rate = %3.2f.',...
                    mupred/mu,pred_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
        else
            update_iter = 1;
        end
    end
    %%
    %%---------------------------------------------------------------
    %% corrector step.
    %%---------------------------------------------------------------
    %%
    if (predcorr) && (~breakyes)
        step_pred = min(pstep,dstep);
        if (mu > 1e-6)
            if (step_pred < 1/sqrt(3));
                expon_used = 1;
            else
                expon_used = max(expon,3*step_pred^2);
            end
        else
            expon_used = max(1,min(expon,3*step_pred^2));
        end
        sigma = min( 1, (mupred/mu)^expon_used );
        sigmu = sigma*mu;
        %%
        if (vers == 1)
            [par,dX,dy,dZ] = HSDHKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
                dX,dZ,coeff,L,X,Z);
        elseif (vers == 2)
            [par,dX,dy,dZ] = HSDNTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
                dX,dZ,coeff,L,X,Z);
        end
        if (solve_ok <= 0)
            msg = 'stop: difficulty in computing corrector directions';
            if (printlevel); fprintf('\n  %s',msg); end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -4;
            break;
        end
        timenew = clock;
        ttime.corr = ttime.corr + etime(timenew,timeold); timeold=timenew;
        %%
        %%-----------------------------------
        %% step-lengths for corrector step
        %%-----------------------------------
        %%
        if (gam == 0)
            gamused = 0.9 + 0.09*min(pstep,dstep);
        else
            gamused = gam;
        end
        kapstep = max( (par.dkap<0)*(-kap/(par.dkap-eps)), (par.dkap>=0)*1e6 );
        taustep = max( (par.dtau<0)*(-tau/(par.dtau-eps)), (par.dtau>=0)*1e6 );
        Xstep = steplength(blk,X,dX,Xchol,invXchol);
        timenew = clock;
        ttime.corr_pstep = ttime.corr_pstep + etime(timenew,timeold); timeold=timenew;
        Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
        timenew = clock;
        pstep = min(1,gamused*min([Xstep,Zstep,kapstep,taustep]));
        dstep = pstep;
        kapcorr = kap + pstep*par.dkap;
        taucorr = tau + pstep*par.dtau;
        trXZcorr = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ)...
            + pstep*dstep*blktrace(blk,dX,dZ);
        mucorr = (trXZcorr+kapcorr*taucorr)/(nn+1);
        ttime.corr_dstep = ttime.corr_dstep + etime(timenew,timeold); timeold=timenew;
        %%
        %%-----------------------------------------
        %%  stopping criteria for corrector step
        %%-----------------------------------------
        %%
        if (iter >= 2)
            idx = max(2,iter-2) : iter;
            corr_slow = all(runhist.gap(idx)./runhist.gap(idx-1) > 0.8);
            idx = max(2,iter-5) : iter;
            corr_convg_rate = mean(runhist.gap(idx)./runhist.gap(idx-1));
            corr_slow = corr_slow + (mucorr/mu > max(min(1,5*corr_convg_rate),0.8));
        end
        if (max(mu,infeas) < 1e-6) && (iter > 10) && (stoplevel) ...
                && (corr_slow && mucorr/mu > 1.0)
            msg = 'stop: lack of progress in corrector';
            if (printlevel)
                fprintf('\n  %s',msg);
                fprintf(': mucorr/mu = %3.2f, corr_convg_rate = %3.2f',...
                    mucorr/mu,corr_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
        else
            update_iter = 1;
        end
    end
    %%
    %%---------------------------------------------------------------
    %% udpate iterate
    %%---------------------------------------------------------------
    %%
    indef = [1 1];
    if (update_iter)
        for t = 1:5
            [Xchol,indef(1)] = blkcholfun(blk,ops(X,'+',dX,pstep));
            timenew = clock;
            ttime.pchol = ttime.pchol + etime(timenew,timeold); timeold = timenew;
            if (indef(1)); pstep = 0.8*pstep; else break; end
        end
        if (t > 1); pstep = gamused*pstep; end
        for t = 1:5
            [Zchol,indef(2)] = blkcholfun(blk,ops(Z,'+',dZ,dstep));
            timenew = clock;
            ttime.dchol = ttime.dchol + etime(timenew,timeold); timeold = timenew;
            if (indef(2)); dstep = 0.8*dstep; else break; end
        end
        if (t > 1); dstep = gamused*dstep; end
        AdX = AXfun(blk,At,par.permA,dX);
        AXtmp = AX(1:m) + pstep*AdX(1:m); tautmp = par.tau+pstep*par.dtau;
        prim_infeasnew = norm(b-AXtmp/tautmp)/normb;
        pinfeas_bad(1) = (prim_infeasnew > max([1e-8,relgap,10*infeas]));
        pinfeas_bad(2) = (prim_infeasnew > max([1e-4,20*prim_infeas]) ...
            && (infeas < 1e-2));
        pinfeas_bad(3) = (max([relgap,dual_infeas]) < 1e-4) ...
            && (prim_infeasnew > max([2*prim_infeas,10*dual_infeas,1e-7]));
        if any(indef)
            msg = 'stop: X, Z not both positive definite';
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -3;
            breakyes = 1;
        elseif any(pinfeas_bad)
            if (stoplevel) && (max(pstep,dstep)<=1) && (kap < 1e-3) ...
                    && (prim_infeasnew > dual_infeas);
                msg = 'stop: primal infeas has deteriorated too much';
                if (printlevel); fprintf('\n  %s, %2.1e',msg,prim_infeasnew);
                    fprintf(' %2.1d,%2.1d,%2.1d',...
                        pinfeas_bad(1),pinfeas_bad(2),pinfeas_bad(3));
                end
                termcode = -7;
                breakyes = 1;
            end
        end
        if (~breakyes)
            X = ops(X,'+',dX,pstep);
            y = y + dstep*dy;  Z = ops(Z,'+',dZ,dstep);
            theta = max(0, theta + pstep*par.dtheta);
            kap = kap + pstep*par.dkap;
            if (tau + pstep*par.dtau > theta); tau = tau + pstep*par.dtau; end
        end
    end
    %%
    %%--------------------------------------------------
    %% perturb Z: do this step before checking for break
    %%--------------------------------------------------
    perturb_Z = 1;
    if (~breakyes) && (perturb_Z)
        trXZtmp = blktrace(blk,X,Z);
        trXE  = blktrace(blk,X,EE);
        Zpert = max(1e-12,0.2*min(relgap,prim_infeas)).*normC./normE;
        Zpert = min(Zpert,0.1*trXZtmp./trXE);
        Zpert = min([1,Zpert,1.5*Zpertold]);
        if (infeas < 1e-2)
            Z = ops(Z,'+',EE,Zpert);
            [Zchol,indef(2)] = blkcholfun(blk,Z);
            if any(indef(2))
                msg = 'stop: Z not positive definite';
                if (printlevel); fprintf('\n  %s',msg); end
                termcode = -3;
                breakyes = 1;
            end
        end
        Zpertold = Zpert;
    end
    %%
    %%---------------------------------------------------------------
    %% compute rp, Rd, infeasibities, etc.
    %%---------------------------------------------------------------
    %%
    y2 = [y; tau; theta];
    AX = AXfun(blk,At,par.permA,X);
    rp = [zeros(m,1); kap; -abar] - AX - Bmat*y2;
    % Rd = ops(Atyfun(blk,At,par.permA,par.isspAy,-y2),'-',Z);
    trXZ = blktrace(blk,X,Z);
    mu  = (trXZ+kap*tau)/(nn+1);
    obj = [blktrace(blk,C,X), b'*y]/tau;
    gap = trXZ/tau^2;
    relgap = gap/(1+mean(abs(obj)));
    ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,[y;0;0]));
    ZpATynorm = ops(ZpATy,'norm');
    prim_infeas = norm(b-AX(1:m)/tau)/normb;
    dual_infeas = ops(ops(C,'-',ops(ZpATy,'/',tau)),'norm')/normC;
    infeas = max(prim_infeas,dual_infeas);
    runhist.pobj(iter+1)    = obj(1);
    runhist.dobj(iter+1)    = obj(2);
    runhist.gap(iter+1)     = gap;
    runhist.relgap(iter+1)  = relgap;
    runhist.pinfeas(iter+1) = prim_infeas;
    runhist.dinfeas(iter+1) = dual_infeas;
    runhist.infeas(iter+1)  = infeas;
    runhist.cputime(iter+1) = etime(clock,tstart);
    runhist.step(iter+1)    = min(pstep,dstep);
    runhist.kappa(iter+1)   = kap;
    runhist.tau(iter+1)     = tau;
    runhist.theta(iter+1)   = theta;
    runhist.useLU(iter+1)   = use_LU;
    timenew = clock;
    ttime.misc = ttime.misc + etime(timenew,timeold); % timeold = timenew;
    [hh,mm,ss] = mytime(sum(runhist.cputime));
    if (printlevel>=3)
        fprintf('\n%2.0f|%4.3f|%4.3f|',iter,pstep,dstep);
        fprintf('%2.1e|%2.1e|%2.1e|',prim_infeas,dual_infeas,gap);
        fprintf('%- 7.6e| %s:%s:%s|',mean(obj),hh,mm,ss);
        fprintf('%2.1e|%2.1e|%2.1e|',kap,tau,theta);
    end
    %%
    %%--------------------------------------------------
    %% check convergence.
    %%--------------------------------------------------
    param.termcode    = termcode;
    param.kap         = kap;
    param.tau         = tau;
    param.theta       = theta;
    param.iter        = iter;
    param.obj         = obj;
    param.gap         = gap;
    param.relgap      = relgap;
    param.mu          = mu;
    param.prim_infeas = prim_infeas;
    param.dual_infeas = dual_infeas;
    param.AX          = AX(1:m)/tau;
    param.ZpATynorm   = ZpATynorm/tau;
    param.normX       = ops(X,'norm')/tau;
    param.normZ       = ops(Z,'norm')/tau;
    if (~breakyes)
        [param,breakyes,use_olditer,msg] = HSDsqlpcheckconvg(param,runhist);
        termcode = param.termcode;  %% important
        if (use_olditer)
            X = ops(X,'-',dX,pstep);
            y = y - dstep*dy;
            Z = ops(Z,'-',dZ,dstep);
            kap   = kap - pstep*par.dkap;
            tau   = tau - pstep*par.dtau;
            theta = theta - pstep*par.dtheta;
            prim_infeas = runhist.pinfeas(iter);
            dual_infeas = runhist.dinfeas(iter);
            gap = runhist.gap(iter); relgap = runhist.relgap(iter);
            obj = [runhist.pobj(iter), runhist.dobj(iter)];
        end
    end
    %%--------------------------------------------------
    %% check for break
    %%--------------------------------------------------
    newtol = max(relgap,infeas);
    update_best(iter+1) = ~( newtol >= besttol ); %#ok
    if update_best(iter+1),
        Xbest = X; ybest = y; Zbest = Z;
        besttol = newtol;
    end
    if besttol < 1e-4 && ~any(update_best(max(1,iter-1):iter+1))
        msg = 'lack of progess in infeas';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -9;
        breakyes = 1;
    end
    if besttol < 1e-3 && newtol > 1.2*besttol && theta < 1e-10 && kap < 1e-6
        msg = 'lack of progress in infeas';
        if (printlevel); fprintf('\n  %s',msg); end
        termcode = -9;
        breakyes = 1;
    end
    if (breakyes > 0.5); break; end
end
%%---------------------------------------------------------------
%% end of main loop
%%---------------------------------------------------------------
%%
use_bestiter = 1;
if (use_bestiter) && (param.termcode <= 0)
    X = Xbest; y = ybest; Z = Zbest;
    % kap = kapbest; tau = taubest; theta = thetabest;
    trXZ = blktrace(blk,X,Z);
    obj  = [blktrace(blk,C,X), b'*y]/tau;
    gap  = trXZ/tau^2;
    relgap = gap/(1+mean(abs(obj)));
    AX     = AXfun(blk,At,par.permA,X);
    ZpATy  = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,[y;0;0]));
    ZpATynorm = ops(ZpATy,'norm');
    prim_infeas = norm(b-AX(1:m)/tau)/normb;
    dual_infeas = ops(ops(C,'-',ops(ZpATy,'/',tau)),'norm')/normC;
    infeas = max(prim_infeas,dual_infeas);
    runhist.pobj(iter+1)    = obj(1);
    runhist.dobj(iter+1)    = obj(2);
    runhist.gap(iter+1)     = gap;
    runhist.relgap(iter+1)  = relgap;
    runhist.pinfeas(iter+1) = prim_infeas;
    runhist.dinfeas(iter+1) = dual_infeas;
    runhist.infeas(iter+1)  = infeas;
end
%%---------------------------------------------------------------
%% produce infeasibility certificates if appropriate
%%---------------------------------------------------------------
%%
X = ops(X,'/',tau); y = y/tau; Z = ops(Z,'/',tau);
if (iter >= 1)
    param.termcode    = termcode;
    param.obj         = obj;
    param.relgap      = relgap;
    param.prim_infeas = prim_infeas;
    param.dual_infeas = dual_infeas;
    param.AX          = AX(1:m)/tau;
    param.ZpATynorm   = ZpATynorm/tau;
    [X,y,Z,resid,reldist,param,msg2] = ...
        HSDsqlpmisc(blk,At,C,b,X,y,Z,par.permZ,param);
    termcode = param.termcode;
end
%%
%%---------------------------------------------------------------
%% recover unrestricted blk from linear blk
%%---------------------------------------------------------------
%%
for p = 1:size(blk,1)
    if (ublkidx(p) == 1)
        n = blk{p,2}/2;
        X{p} = X{p}(1:n)-X{p}(n+1:2*n);
        Z{p} = Z{p}(1:n);
    end
end
%%
%%---------------------------------------------------------------
%%  print summary
%%---------------------------------------------------------------
%%
dimacs = [prim_infeas; 0; dual_infeas; 0];
dimacs = [dimacs; [-diff(obj); gap]/(1+sum(abs(obj)))];
info.dimacs   = dimacs;
info.termcode = termcode;
info.iter     = iter;
info.obj      = obj;
info.gap      = gap;
info.relgap   = relgap;
info.pinfeas  = prim_infeas;
info.dinfeas  = dual_infeas;
info.cputime  = sum(runhist.cputime);
info.time     = ttime;
info.resid    = resid;
info.reldist  = reldist;
info.normX    = ops(X,'norm');
info.normy    = norm(y);
info.normZ    = ops(Z,'norm');
info.normA    = ops(At,'norm');
info.normb    = norm(b);
info.normC    = ops(C,'norm');
info.msg1     = msg;
info.msg2     = msg2;
%%
if isoctave,
  warning(w1.state,w1.identifier);
else
  warning(w2.state,w2.identifier);
  warning(w1.state,w1.identifier);
end
sqlpsummary(info,ttime,[],printlevel);
%%*****************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function flip = spm_flip_analyze_images
% Do Analyze format images need to be left-right flipped?
%-----------------------------------------------------------------------
% @(#)spm_flip_analyze_images.m	2.4 03/05/07

global defaults
if isempty(defaults) | ~isfield(defaults,'analyze') |...
         ~isfield(defaults.analyze,'flip')
	warning(['Cant get default Analyze orientation - assuming' ...
		 ' not flipped']);
	flip = 0;
	return;
end;
flip = defaults.analyze.flip;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%   4. Mass Points
%   5. Porous Points
%	6. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%   7. 3-Element Hill Muscle Model
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting-lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Actual DRIVER of the code, where the time-stepping occurs ->
%           gets called by main2d to do the "magic" :)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X, Y, U, V, xLag, yLag] = IBM_Driver(struct_name, mu, rho, grid_Info, dt, T_FINAL, model_Info)

%
%    2D IMMERSED BOUNDARY SOLVER ON RECTANGULAR DOMAIN w/ PERIODIC BOUNDARIES
%
%    x-Momentum Conservation: rho*u_t = -rho*u*u_x + rho*v*u_y + mu*laplacian(u) - p_x + F_x
%    y-Momentum Convervation: rho*v_t = -rho*u*v_x + rho*v*v_y + mu*laplacian(v) - p_y + F_y
%
%    Incompressibility: u_x + v_y = 0
%
%    LagPts/dt = int{ u(x,t) delta( x - LagPts(s,t) ) dx }
%    F_x = int{ fx(s,t) delta(x - LagPts(s,t)) ds }
%    F_y = int{ fy(s,t) delta(x - LagPts(s,t)) ds }
%
fprintf('\n________________________________________________________________________________\n\n');
fprintf('\n---------------->>                 IB2d                      <<----------------\n');
fprintf('\n________________________________________________________________________________\n\n');
fprintf('If using the code for research purposes please cite the following two papers: \n');
fprintf('     [1] N.A. Battista, A.J. Baird, L.A. Miller, A mathematical model and MATLAB code for muscle-fluid-structure simulations, Integ. Comp. Biol. 55(5):901-11 (2015)\n');
fprintf('     [2] N.A. Battista, W.C. Strickland, L.A. Miller, IB2d a Python and MATLAB implementation of the immersed boundary method, Bioinspir. Biomim. 12(3):036003 (2017)');
fprintf('\n________________________________________________________________________________');

fprintf('\n\n\n |****** Prepping Immersed Boundary Simulation ******|\n');
fprintf('\n\n--> Reading input data for simulation...\n\n');

% Temporal Information
NTime = floor(T_FINAL/dt)+1; % # of total time-steps (floor'd so exact number of time-steps)
dt = T_FINAL/NTime;          %time-step (slightly perturbed dt, so exact # of time-steps are used)
current_time = 0.0;          %initialize start of simulation to time, 0

% GRID INFO %
Nx =   grid_Info(1); % # of Eulerian pts. in x-direction
Ny =   grid_Info(2); % # of Eulerian pts. in y-direction
Lx =   grid_Info(3); % Length of Eulerian grid in x-coordinate
Ly =   grid_Info(4); % Length of Eulerian grid in y-coordinate
dx =   grid_Info(5); % Spatial-size in x
dy =   grid_Info(6); % Spatial-size in y
supp = grid_Info(7); % Delta-function support

% PRINTING/PLOTTING INFO %
pDump = grid_Info(8);           % Print (Plot) Dump interval
pMatlab = grid_Info(9);         % Plot in Matlab? (1=YES,0=NO)
lagPlot = grid_Info(10);        % Plot LAGRANGIAN PTs ONLY in Matlab
velPlot = grid_Info(11);        % Plot LAGRANGIAN PTs + VELOCITY FIELD in Matlab
vortPlot = grid_Info(12);       % Plot LAGRANGIAN PTs + VORTICITY colormap in Matlab
uMagPlot = grid_Info(13);       % Plot LAGRANGIAN PTs + MAGNITUDE OF VELOCITY colormap in Matlab
pressPlot = grid_Info(14);      % Plot LAGRANGIAN PTs + PRESSURE colormap in Matlab

% MODEL STRUCTURE DATA STORED %
springs_Yes = model_Info(1);           % Springs: 0 (for no) or 1 (for yes)
update_Springs_Flag = model_Info(2);   % Update_Springs: 0 (for no) or 1 (for yes)
target_pts_Yes = model_Info(3);        % Target_Pts: 0 (for no) or 1 (for yes)
update_Target_Pts = model_Info(4);     % Update_Target_Pts: 0 (for no) or 1 (for yes)
beams_Yes = model_Info(5);             % Beams: 0 (for no) or 1 (for yes)
update_Beams_Flag = model_Info(6);     % Update_Beams: 0 (for no) or 1 (for yes)
nonInv_beams_Yes = model_Info(7);      % Beams (non-invariant): 0 (for no) or 1 (for yes)
update_nonInv_Beams_Flag = model_Info(8); % Update_nonInv_Beams: 0 (for no) or 1 (for yes)
muscles_Yes = model_Info(9);           % FV-LT Muscles: 0 (for no) or 1 (for yes)
hill_3_muscles_Yes = model_Info(10);    % Hill 3-Element Muscle: 0 (for no) or 1 (for yes)
arb_ext_force_Yes = model_Info(11);     % Arbitrary External Force: 0 (for no) or 1 (for yes)
tracers_Yes = model_Info(12);          % Tracers: 0 (for no) or 1 (for yes)
mass_Yes = model_Info(13);             % Mass Points: 0 (for no) or 1 (for yes)
gravity_Yes = model_Info(14);          % Gravity: 0 (for no) or 1 (for yes)
%NOTE: model_Info(15)/(16) - components of gravity vector
porous_Yes = model_Info(17);           % Porous Media: 0 (for no) or 1 (for yes)
concentration_Yes = model_Info(18);    % Background Concentration Gradient: 0 (for no) or 1 (for yes)
electro_phys_Yes = model_Info(19);     % Electrophysiology (FitzHugh-Nagumo): 0 (for no) or 1 (for yes)
d_Springs_Yes = model_Info(20);        % Damped Springs: 0 (for no) or 1 (for yes)
update_D_Springs_Flag = model_Info(21);% Update_Damped_Springs: % 0 (for no) or 1 (for yes)
boussinesq_Yes = model_Info(22);       % Boussinesq Approx.: 0 (for no) or 1 (for yes)
exp_Coeff = model_Info(23);            % Expansion Coefficient (e.g., thermal, etc) for Boussinesq approx.
general_force_Yes = model_Info(24);    % General User-Defined Force Term: 0 (for no) or 1 (for yes)

%Lagrangian Structure Data
ds = Lx / (2*Nx);                   %Lagrangian Spacing
grid_Info(9) = ds;

% Create EULERIAN Mesh (these assume periodicity in x and y)
x = (0:dx:Lx-dx);  X = [];
y = (0:dy:Ly-dy)'; Y = [];
%Create x-Mesh
for i=1:Nx
    X = [X; x];
end
%Create y-Mesh
for i=1:Ny
    Y = [Y y];
end

% % % % % HOPEFULLY WHERE I CAN READ IN INFO!!! % % % % %

% READ IN LAGRANGIAN POINTS %
[Nb,xLag,yLag] = read_Vertex_Points(struct_name);
grid_Info(8) = Nb;          % # Total Number of Lagrangian Pts.
xLag_P = xLag;              % Initialize previous Lagrangian x-Values (for use in muscle-model)
yLag_P = yLag;              % Initialize previous Lagrangian y-Values (for use in muscle-model)

fprintf('\n--> FIBER MODEL INCLUDES: \n');

% READ IN SPRINGS (IF THERE ARE SPRINGS) %
if ( springs_Yes == 1 )
    fprintf('  -Springs and ');
    if update_Springs_Flag == 0
        fprintf('NOT dynamically updating spring properties\n');
    else
        fprintf('dynamically updating spring properties\n');
    end
    springs_info = read_Spring_Points(struct_name);
        %springs_info: col 1: starting spring pt (by lag. discretization)
        %              col 2: ending spring pt. (by lag. discretization)
        %              col 3: spring stiffness
        %              col 4: spring resting lengths
        %              col 5: spring linearity (1=linear, >1 non-linear)
else
    springs_info = 0;  %just to pass placeholder into "please_Find_Lagrangian_Forces_On_Eulerian_grid function"
end

% READ IN BEAMS (IF THERE ARE BEAMS aka TORSIONAL SPRINGS) %
if ( beams_Yes == 1)
    fprintf('  -Beams ("Torsional Springs") and ');
    if update_Beams_Flag == 0
        fprintf('NOT dynamically updating beam properties\n');
    else
        fprintf('dynamically updating beam properties\n');
    end

    beams_info = read_Beam_Points(struct_name);
    %beams:      col 1: 1ST PT.
    %            col 2: MIDDLE PT. (where force is exerted)
    %            col 3: 3RD PT.
    %            col 4: beam stiffness
    %            col 5: curavture
else
    beams_info = 0;
end

% READ IN BEAMS (IF THERE ARE NONINVARIANT BEAMS) %
if ( nonInv_beams_Yes == 1)
    fprintf('  -Beams ("non-Invariant") and ');
    if update_nonInv_Beams_Flag == 0
        fprintf('NOT dynamically updating non-invariant beam properties\n');
    else
        fprintf('dynamically updating non-invariant beam properties\n');
    end

    nonInv_beams_info = read_nonInv_Beam_Points(struct_name);
    %beams:      col 1: 1ST PT.
    %            col 2: MIDDLE PT. (where force is exerted)
    %            col 3: 3RD PT.
    %            col 4: beam stiffness
    %            col 5: x-curavture
    %            col 6: y-curavture
else
    nonInv_beams_info = 0;
end

% READ IN TARGET POINTS (IF THERE ARE TARGET PTS) %
if ( target_pts_Yes == 1)
    fprintf('  -Target Pts. and ');
    if update_Target_Pts == 0
        fprintf('NOT dynamically updating target point properties\n');
    else
        fprintf('dynamically updating target point properties\n');
    end

    target_aux = read_Target_Points(struct_name);
    %target_aux: col 1: Lag Pt. ID w/ Associated Target Pt.
    %            col 2: target STIFFNESSES
    target_info(:,1) = target_aux(:,1); %Stores Lag-Pt IDs in col vector
    for i=1:length(target_info(:,1))
        id = target_info(i,1);
        target_info(i,2) = xLag(id);    %Stores Original x-Lags as x-Target Pt. Identities
        target_info(i,3) = yLag(id);    %Stores Original y-Lags as y-Target Pt. Identities
    end

    target_info(:,4) = target_aux(:,2); %Stores Target Stiffnesses
else
    target_info = 0;
end

% READ IN MASS POINTS (IF THERE ARE MASS PTS) %
if ( mass_Yes == 1)
    fprintf('  -Mass Pts. with ');
    if gravity_Yes == 0
        fprintf('NO artificial gravity\n');
    else
        fprintf('artificial gravity\n');
    end
    mass_aux = read_Mass_Points(struct_name);
    %target_aux: col 1: Lag Pt. ID w/ Associated Mass Pt.
    %            col 2: "Mass-spring" stiffness parameter
    %            col 3: "MASS" value parameter
    mass_info(:,1) = mass_aux(:,1); %Stores Lag-Pt IDs in col vector
    for i=1:length(mass_info(:,1))
        id = mass_info(i,1);
        mass_info(i,2) = xLag(id);    %Stores Original x-Lags as x-Mass Pt. Identities
        mass_info(i,3) = yLag(id);    %Stores Original y-Lags as y-Mass Pt. Identities
    end

    mass_info(:,4) = mass_aux(:,2);   %Stores "mass-spring" parameter
    mass_info(:,5) = mass_aux(:,3);   %Stores "MASS" value parameter
else
    mass_info = 0;
end

% CONSTRUCT GRAVITY INFORMATION (IF THERE IS GRAVITY) %
if gravity_Yes == 1
    %gravity_Vec(1) = model_Info(12);     % x-Component of Gravity Vector
    %gravity_Vec(2) = model_Info(13);     % y-Component of Gravity Vector
    xG = model_Info(15);
    yG = model_Info(16);
    normG = sqrt( xG^2 + yG^2 );
    gravity_Info = [gravity_Yes xG/normG yG/normG];
    %   col 1: flag if considering gravity
    %   col 2: x-component of gravity vector (normalized)
    %   col 3: y-component of gravity vector (normalized)

    clear xG yG normG;

else
    gravity_Info = 0;
end

% READ IN POROUS MEDIA INFO (IF THERE IS POROSITY) %
if ( porous_Yes == 1)
    fprintf('  -Porous Points\n');
    porous_aux = read_Porous_Points(struct_name);
    %porous_aux: col 1: Lag Pt. ID w/ Associated Porous Pt.
    %            col 2: Porosity coefficient
    %            col 3: Flag for derivative stencil!
    porous_info(:,1) = porous_aux(:,1); %Stores Lag-Pt IDs in col vector
    for i=1:length(porous_info(:,1))
        id = porous_info(i,1);
        porous_info(i,2) = xLag(id);    %Stores Original x-Lags as x-Porous Pt. Identities
        porous_info(i,3) = yLag(id);    %Stores Original y-Lags as y-Porous Pt. Identities
    end

    porous_info(:,4) = porous_aux(:,2); %Stores Porosity Coefficient
    porous_info(:,5) = porous_aux(:,3); %Stores flag for derivative stencil
else
    porous_info = 0;
end

% READ IN DAMPED SPRINGS (IF THERE ARE SPRINGS) %
if ( d_Springs_Yes == 1 )
    fprintf('  -Damped Springs and ');
    if update_D_Springs_Flag == 0
        fprintf('NOT dynamically updating damped spring properties\n');
    else
        fprintf('dynamically updating spring properties\n');
    end
    d_springs_info = read_Damped_Spring_Points(struct_name);
        %springs_info: col 1: starting spring pt (by lag. discretization)
        %              col 2: ending spring pt. (by lag. discretization)
        %              col 3: spring stiffness
        %              col 4: spring resting lengths
        %              col 5: damping coefficients
else
    d_springs_info = 0;  %just to pass placeholder into "please_Find_Lagrangian_Forces_On_Eulerian_grid function"
end

% READ IN MUSCLES (IF THERE ARE MUSCLES) %
if ( muscles_Yes == 1 )
    fprintf('  -MUSCLE MODEL (Force-Velocity / Length-Tension Model)\n');
    muscles_info = read_Muscle_Points(struct_name);
        %         muscles: col 1: MASTER NODE (by lag. discretization)
        %         col 2: SLAVE NODE (by lag. discretization)
        %         col 3: length for max. muscle tension
        %         col 4: muscle constant
        %         col 5: hill parameter, a
        %         col 6: hill parameters, b
        %         col 7: force maximum!
else
    muscles_info = 0;  %just to pass placeholder into "please_Find_Lagrangian_Forces_On_Eulerian_grid function"
end

% READ IN MUSCLES (IF THERE ARE MUSCLES) %
if ( hill_3_muscles_Yes == 1 )
    fprintf('  -MUSCLE MODEL (3 Element Hill Model)\n');
    muscles3_info = read_Hill_3Muscle_Points(struct_name);
        %         muscles: col 1: MASTER NODE (by lag. discretization)
        %         col 2: SLAVE NODE (by lag. discretization)
        %         col 3: length for max. muscle tension
        %         col 4: muscle constant
        %         col 5: hill parameter, a
        %         col 6: hill parameters, b
        %         col 7: force maximum!
        %         col 8: kStiffness of NL Spring
        %         col 9: alpha, degree of non-linearity
else
    muscles3_info = 0;  %just to pass placeholder into "please_Find_Lagrangian_Forces_On_Eulerian_grid function"
end

% READ IN GENERAL FORCE PARAMETERS (IF THERE IS A USER-DEFINED FORCE) %
if ( general_force_Yes == 1 )
    fprintf('  -GENERAL FORCE MODEL (user-defined force term)\n');
    gen_force_info = read_General_Forcing_Function(struct_name);
        %
        %
        % ALL PARAMETERS / FORCE FUNCTION SET BY USER!
        %
        %
else
    gen_force_info = 0;  %just to pass placeholder into "please_Find_Lagrangian_Forces_On_Eulerian_grid function"
end

% SOLVE ELECTROPHYSIOLOGY MODEL FOR PUMPING %
if electro_phys_Yes == 1
    fprintf('  -Electrophysiology model via FitzHugh-Nagumo\n');
    fprintf('\n\n--> Solving Electrophysiology Model...\n');
    [electro_potential, ePhys_Start, ePhys_End] = FitzHugh_Nagumo_1d(dt,T_FINAL);
    ePhys_Ct = 1;
    fprintf('--> Finished Computing Electrophysiology...time for IBM!\n');
else
    electro_potential = 0;
end

fprintf('\n\n--> Background Flow Items\n');
if ( tracers_Yes == 0 ) && (concentration_Yes == 0)
    fprintf('      (No tracers nor other passive scalars immersed in fluid)\n\n');
end

% READ IN TRACERS (IF THERE ARE TRACERS) %
if (tracers_Yes == 1)
   fprintf('  -Tracer Particles included\n');
   [~,xT,yT] = read_Tracer_Points(struct_name);
   tracers = zeros(length(xT),4);
   tracers(1,1) = 1;
   tracers(:,2) = xT;
   tracers(:,3) = yT;
        %tracers_info: col 1: xPt of Tracers
        %              col 2: yPt of Tracers
else
   tracers = 0;
end

% READ IN CONCENTRATION (IF THERE IS A BACKGROUND CONCENTRATION) %
if ( concentration_Yes == 1 )
    fprintf('  -Background concentration included\n');
    [C,kDiffusion] = read_In_Concentration_Info(struct_name,Nx);
        %C:           Initial background concentration
        %kDiffusion:  Diffusion constant for Advection-Diffusion
else
    C = 0; % placeholder for plotting
end

% CONSTRUCT BOUSSINESQ INFORMATION (IF USING BOUSSINESQ) %
if boussinesq_Yes == 1

    if length(gravity_Info) == 1
        fprintf('\n\n\n READ THE ERROR MESSAGE -> YOU MUST FLAG GRAVITY IN INPUT FILE FOR BOUSSINESQ! :) \n\n\n');
        error('YOU MUST FLAG GRAVITY IN INPUT FILE FOR BOUSSINESQ! :)');
    elseif concentration_Yes == 0
        fprintf('\n\n\n READ THE ERROR MESSAGE -> YOU MUST HAVE BACKGROUND CONCENTRATION FOR BOUSSINESQ! :) \n\n\n');
        error('YOU MUST FLAG CONCENTRATION IN INPUT FILE FOR BOUSSINESQ! :)');
    end

    % Forms Boussinesq forcing terms, e.g., (exp_Coeff)*gVector for Momentum Eq.
    [fBouss_X,fBouss_Y] = please_Form_Boussinesq_Forcing_Terms(exp_Coeff,Nx,Ny,gravity_Info);

    % Finds initial concentration Laplacian
    Cxx = DD(C,dx,'x');
    Cyy = DD(C,dy,'y');
    laplacian_C = Cxx+Cyy;
    C_0 = zeros(size(C)); % Define background concentration
end

% Initialize the initial velocities to zero.
U = zeros(Ny,Nx);                                % x-Eulerian grid velocity
V = U;                                           % y-Eulerian grid velocity
mVelocity = zeros( length(mass_info(:,1)), 2 );  % mass-Pt velocity

if arb_ext_force_Yes == 1
    firstExtForce = 1;                           % initialize external forcing
    indsExtForce = 0;                            % initialize for external forcing computation
end

% ACTUAL TIME-STEPPING IBM SCHEME!
%(flags for storing structure connects for printing and printing to .vtk)
cter = 0; ctsave = 0; firstPrint = 1; loc = 1; diffy = 1;

% CREATE VIZ_IB2D FOLDER and VISIT FILES
mkdir('viz_IB2d');
mkdir('hier_IB2d_data');
%cd('viz_IB2d');
%vizID = fopen('dumps.visit','w');
%fprintf(vizID,'!NBLOCKS 6\n');
%cd ..

%Initialize Vorticity, uMagnitude, and Pressure for initial colormap
%Print initializations to .vtk
vort=zeros(Nx,Ny); uMag=vort; p = vort;  lagPts = [xLag yLag zeros(length(xLag),1)];
[connectsMat,spacing] = give_Me_Lag_Pt_Connects(ds,xLag,yLag,Nx,springs_Yes,springs_info);
Fxh = vort; Fyh = vort; F_Lag = zeros( length(xLag), 2);
print_vtk_files(ctsave,vort,uMag,p,U,V,Lx,Ly,Nx,Ny,lagPts,springs_Yes,connectsMat,tracers,concentration_Yes,C,Fxh,Fyh,F_Lag);
fprintf('\n |****** Begin IMMERSED BOUNDARY SIMULATION! ******| \n\n');
fprintf('Current Time(s): %6.6f\n\n',current_time);
ctsave = ctsave+1;

%
%
%
% **************************************************************************
% * * * * * * * * * * * * * BEGIN TIME-STEPPING! * * * * * * * * * * * * * *
% **************************************************************************
%
%
%
while current_time < T_FINAL

    %
    %
    %**************** Step 1: Update Position of Boundary of membrane at half time-step *******************
    %                           (Variables end with h if it is a half-step)
    %
    [xLag_h, yLag_h] = please_Move_Lagrangian_Point_Positions(U, V, xLag, yLag, xLag, yLag, x, y, dt/2, grid_Info,0);

    if mass_Yes == 1
       [mass_info, massLagsOld] = please_Move_Massive_Boundary(dt/2,mass_info,mVelocity);
    end

    if ( ( electro_phys_Yes == 1) && (muscles_Yes == 0) )
        springs_info(ePhys_Start:ePhys_End,3) = ( 8.5e0*electro_potential(ePhys_Ct,:)') .^4;%( 1e4*electro_potential(ePhys_Ct,:)'.*springs_info(ePhys_Start:ePhys_End,3) ).^4;
        ePhys_Ct = ePhys_Ct + 1;
    end

    if ( ( update_Springs_Flag == 1 ) && ( springs_Yes == 1 ) )
       springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info);
    end

    if ( ( update_Target_Pts == 1 ) && ( target_pts_Yes == 1) )
       target_info = update_Target_Point_Positions(dt,current_time,target_info);
    end

    if ( ( update_Beams_Flag == 1 ) && ( beams_Yes == 1) )
       beams_info = update_Beams(dt,current_time,beams_info);
    end

    if ( ( update_nonInv_Beams_Flag == 1 ) && ( nonInv_beams_Yes == 1) )
       nonInv_beams_info = update_nonInv_Beams(dt,current_time,nonInv_beams_info);
    end

    if ( ( update_D_Springs_Flag == 1 ) && ( d_Springs_Yes == 1) )
       d_springs_info = update_Damped_Springs(dt,current_time,d_springs_info);
    end

    %
    %
    %**************** STEP 2: Calculate Force coming from membrane at half time-step ****************
    %
    %
    [Fxh, Fyh, F_Mass_Bnd, F_Lag] =    please_Find_Lagrangian_Forces_On_Eulerian_grid(dt, current_time, xLag_h, yLag_h, xLag_P, yLag_P, x, y, grid_Info, model_Info, springs_info, target_info, beams_info, nonInv_beams_info ,muscles_info, muscles3_info, mass_info, electro_potential, d_springs_info, gen_force_info);

    % Once force is calculated, can finish time-step for massive boundary
    if mass_Yes == 1
        % Update Massive Boundary Velocity
        mVelocity_h = please_Update_Massive_Boundary_Velocity(dt/2,mass_info,mVelocity,F_Mass_Bnd,gravity_Info);

        % Update Massive Boundary Position for Time-step
        mass_info(:,[2 3]) = massLagsOld;
        [mass_info,~] = please_Move_Massive_Boundary(dt,mass_info,mVelocity_h);

        % Update Massive Boundary Velocity for Time-step
        mVelocity = please_Update_Massive_Boundary_Velocity(dt,mass_info,mVelocity,F_Mass_Bnd,gravity_Info);
    end

    % Add artificial force from fluid boundary, if desired.
    if arb_ext_force_Yes == 1
        [Fx_Arb, Fy_Arb, firstExtForce, indsExtForce] = please_Compute_External_Forcing(dt, current_time, x, y, grid_Info, U, V, firstExtForce, indsExtForce);
        Fxh = Fxh + Fx_Arb;
        Fyh = Fyh + Fy_Arb;
    end

    %
    %
    %***************** STEP 3: Solve for Fluid motion ******************************************
    %
    %
    % Add in effect from BOUSSINESQ
    if boussinesq_Yes == 1
        Fxh = Fxh + rho*fBouss_X*(C);
        Fyh = Fyh + rho*fBouss_Y*(C);
        [Uh, Vh, U, V, p] =   please_Update_Fluid_Velocity(U, V, Fxh, Fyh, rho, mu, grid_Info, dt);
    else
        [Uh, Vh, U, V, p] =   please_Update_Fluid_Velocity(U, V, Fxh, Fyh, rho, mu, grid_Info, dt);
    end

    %
    %
    %***************** STEP 4: Update Position of Boundary of membrane again for a half time-step *******
    %
    %

    xLag_P = xLag_h;   % Stores old Lagrangian x-Values (for muscle model)
    yLag_P = yLag_h;   % Stores old Lagrangian y-Values (for muscle model)
    %Uh, Vh instead of U,V?
    [xLag, yLag] =     please_Move_Lagrangian_Point_Positions(Uh, Vh, xLag, yLag, xLag_h, yLag_h, x, y, dt, grid_Info,porous_Yes);

    %NOTE: SET UP FOR BOTH CLOSED + OPEN SYSTEMS NOW!!!
    if porous_Yes == 1
       [Por_Mat,nX,nY] = please_Compute_Porous_Slip_Velocity(ds,xLag,yLag,porous_info,F_Lag);
       xLag( porous_info(:,1) ) = xLag( porous_info(:,1) ) - dt*( Por_Mat(:,1)+Por_Mat(:,2) ).*nX;
       yLag( porous_info(:,1) ) = yLag( porous_info(:,1) ) - dt*( Por_Mat(:,1)+Por_Mat(:,2) ).*nY;
       porous_info(:,2) = xLag( porous_info(:,1) );    %Stores x-Lags as x-Porous Pt. Identities
       porous_info(:,3) = yLag( porous_info(:,1) );    %Stores y-Lags as y-Porous Pt. Identities
       xLag = mod(xLag, Lx); % If structure goes outside domain
       yLag = mod(yLag, Ly); % If structure goes outside domain
    end

    % If there are tracers, update tracer positions %
    if tracers_Yes == 1
        %Uh, Vh instead of U,V?
        [xT, yT] = please_Move_Lagrangian_Point_Positions(Uh, Vh, xT, yT, xT, yT, x, y, dt, grid_Info,0); %0 for always no porous tracers
        tracers(:,2) = xT;
        tracers(:,3) = yT;
    end

    % If there is a background concentration, update concentration-gradient %
    if concentration_Yes == 1
       %[C,~] = please_Update_Adv_Diff_Concentration_Flux_Limiter_FV(C,dt,dx,dy,U,V,kDiffusion);
       [C,~] = please_Update_Adv_Diff_Concentration(C,dt,dx,dy,U,V,kDiffusion);
    end

    % Plot Lagrangian/Eulerian Dynamics! %
    if ( ( mod(cter,pDump) == 0  ) && ( cter >= pDump ) )

        %Compute vorticity, uMagnitude
        vort = give_Me_Vorticity(U,V,dx,dy);
        uMag = give_Me_Magnitude_Velocity(U,V);

        %Plot in Matlab
        if pMatlab == 1
            [loc, diffy] = please_Plot_Results(ds,X,Y,U,V,vort,uMag,p,xLag,yLag,lagPlot,velPlot,vortPlot,pressPlot,uMagPlot,firstPrint,loc,diffy,spacing);
        end

        %Print .vtk files!
        lagPts = [xLag yLag zeros(length(xLag),1)];
        print_vtk_files(ctsave,vort,uMag',p',U',V',Lx,Ly,Nx,Ny,lagPts,springs_Yes,connectsMat,tracers,concentration_Yes,C,Fxh',Fyh',F_Lag);

        %Print Current Time
        fprintf('Current Time(s): %6.6f\n',current_time);

        % Prints CFL For Advection:
        maxy = max( max(max(abs(U))), max(max(abs(V))) );
        CFL_adv = (dt/dx)*maxy;
        %CFL_para = dt/dx^2*k;
        %fprintf('CFL: %d\n\n',max(CFL_adv,CFL_para));
        fprintf('CFL: %d\n\n',CFL_adv);

        %Update print counter for filename index
        ctsave=ctsave+1; firstPrint = 0;
        %pause();
         %electro_potential(ePhys_Ct,:)'
         %springs_info(ePhys_Start:ePhys_End,3)
         %springs_info(:,3)
         %pause();

    end

    % Update current_time value & counter
    current_time = current_time+dt;
    cter = cter + 1;
    %pause(0.);

end %ENDS TIME-STEPPING LOOP

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives appropriate string number for filename in printing the
% .vtk files. for viz_IB2d and hier_IB2d_data folders
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_vtk_files(ctsave,vort,uMag,p,U,V,Lx,Ly,Nx,Ny,lagPts,springs_Yes,connectsMat,tracers,concentration_Yes,C,fXGrid,fYGrid,F_Lag)

%Give spacing for grid
dx = Lx/Nx;
dy = Ly/Ny;

%
% Print all Eulerian Data and Lagrangian Position Data to viz_IB2d folder
%
cd('viz_IB2d'); %Go into viz_IB2d directory

    %Find string number for storing files
    strNUM = give_String_Number_For_VTK(ctsave);
    vortfName = ['Omega.' strNUM '.vtk'];
    uMagfName = ['uMag.' strNUM '.vtk'];
    pfName = ['P.' strNUM '.vtk'];
    uXName = ['uX.' strNUM '.vtk'];
    uYName = ['uY.' strNUM '.vtk'];
    fXName = ['fX.' strNUM '.vtk'];
    fYName = ['fY.' strNUM '.vtk'];
    fMagName = ['fMag.' strNUM '.vtk'];
    velocityName = ['u.' strNUM '.vtk'];
    lagPtsName = ['lagsPts.' strNUM '.vtk'];

    %Print Lagrangian Pts to .vtk format
    savevtk_points(lagPts, lagPtsName, 'lagPts');

    % Print Spring Connections (* if springs *)
    if springs_Yes == 1
         %Print Lagrangian Pts w/ CONNECTIONS to .vtk format
        lagPtsConName=['lagPtsConnect.' strNUM '.vtk'];
        savevtk_points_connects(lagPts, lagPtsConName, 'lagPtsConnected',connectsMat);
    end

    %Print Tracer Pts (*if tracers*)
    if tracers(1,1) == 1
        tracersPtsName = ['tracer.' strNUM '.vtk'];
        %tMatrix = tracers(:,2:4);
        savevtk_points(tracers(:,2:4),tracersPtsName, 'tracers');
    end

    %Print another cycle to .visit file
    %fprintf(vizID,[vortfName '\n']);
    %fprintf(vizID,[uMagfName '\n']);
    %fprintf(vizID,[pfName '\n']);
    %fprintf(vizID,[uXName '\n']);
    %fprintf(vizID,[uYName '\n']);
    %fprintf(vizID,[velocityName '\n']);

    %Print SCALAR DATA (i.e., colormap data) to .vtk file
    savevtk_scalar(vort, vortfName, 'Omega',dx,dy);
    savevtk_scalar(uMag, uMagfName, 'uMag',dx,dy);
    savevtk_scalar(p, pfName, 'P',dx,dy);
    savevtk_scalar(U, uXName, 'uX',dx,dy);
    savevtk_scalar(V, uYName, 'uY',dx,dy);
    savevtk_scalar(fXGrid, fXName, 'fX',dx,dy);
    savevtk_scalar(fYGrid, fYName, 'fY',dx,dy);
    savevtk_scalar(sqrt( fXGrid.^2 + fYGrid.^2 ), fMagName, 'fMag',dx,dy);

    if concentration_Yes == 1
        confName = ['concentration.' strNUM '.vtk'];
        savevtk_scalar(C', confName, 'Concentration',dx,dy);
    end

    %Print VECTOR DATA (i.e., velocity data) to .vtk file
    savevtk_vector(U, V, velocityName, 'u',dx,dy)

%Get out of viz_IB2d folder
cd ..

if length( lagPts ) <= 5
    %
    % Print Lagrangian Force Data to hier_IB2d_data folder (if <= 5 lag pts)
    %
    cd('hier_IB2d_data'); %change directory to hier-data folder

        % Save x-y force data!
        fLag_XName = ['fX_Lag.' strNUM '.vtk'];
        fLag_YName = ['fY_Lag.' strNUM '.vtk'];
        savevtk_points_with_scalar_data( lagPts, F_Lag(:,1), fLag_XName, 'fX_Lag');
        savevtk_points_with_scalar_data( lagPts, F_Lag(:,2), fLag_YName, 'fY_Lag');

        % Save force Magnitude (no normal/tangential -> not enough points)
        fMagName = ['fMag.' strNUM '.vtk'];
        fLagMag = sqrt( F_Lag(:,1).^2 + F_Lag(:,2).^2 ); % Compute magnitude of forces on boundary
        savevtk_points_with_scalar_data( lagPts, fLagMag, fMagName, 'fMag');
    cd ..
else

    %
    % Print Lagrangian Force Data to hier_IB2d_data folder
    %
    [F_Tan_Mag,F_Normal_Mag] = please_Compute_Normal_Tangential_Forces_On_Lag_Pts(lagPts,F_Lag);
    %
    cd('hier_IB2d_data'); %change directory to hier-data folder

        % Save x-y force data!
        fLag_XName = ['fX_Lag.' strNUM '.vtk'];
        fLag_YName = ['fY_Lag.' strNUM '.vtk'];
        savevtk_points_with_scalar_data( lagPts, F_Lag(:,1), fLag_XName, 'fX_Lag');
        savevtk_points_with_scalar_data( lagPts, F_Lag(:,2), fLag_YName, 'fY_Lag');

        % Save force Magnitude, Mag. Normal, Mag. Tangential
        fMagName = ['fMag.' strNUM '.vtk'];
        fNormalName = ['fNorm.' strNUM '.vtk'];
        fTangentName = ['fTan.' strNUM '.vtk'];

        fLagMag = sqrt( F_Lag(:,1).^2 + F_Lag(:,2).^2 ); % Compute magnitude of forces on boundary

        savevtk_points_with_scalar_data( lagPts, fLagMag, fMagName, 'fMag');
        savevtk_points_with_scalar_data( lagPts, F_Normal_Mag, fNormalName, 'fNorm');
        savevtk_points_with_scalar_data( lagPts, F_Tan_Mag, fTangentName, 'fTan');

    cd .. % Get out of hier_IB2d_data folder

end % ENDS IF-STATEMENT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: give me Connects Vector for printing Lagrangian .vtk info!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [connectsMat,space] = give_Me_Lag_Pt_Connects(ds,xLag,yLag,Nx,springs_Yes,springs_info)

        %springs_info: col 1: starting spring pt (by lag. discretization)
        %              col 2: ending spring pt. (by lag. discretization)
        %              col 3: spring stiffness
        %              col 4: spring resting lengths
        %              col 5: spring linearity (1=linear, >1 non-linear)

space = 20*ds;

if springs_Yes == 1
    connectsMat(:,1) = springs_info(:,1)-1; % -1 For Cpp notation (and .vtk counting)
    connectsMat(:,2) = springs_info(:,2)-1; % -1 For Cpp notation (and .vtk counting)
else

    connectsMat = 0;
end

%N = length(xLag);

% if Nx <= 32
%     space = 5*ds;
% elseif Nx <= 64
%    space = 5*ds;
% elseif Nx <=128
%    space = 5*ds;
% elseif Nx <=256
%     space = 10*ds;
% elseif Nx <= 512
%     space = 20*ds;
% else
%     space = 40*ds;
% end

% ct = 1;
% for i=1:N
%     if i<N
%         x1=xLag(i); x2=xLag(i+1);
%         y1=yLag(i); y2=yLag(i+1);
%         dist = sqrt( (x1-x2)^2 + (y1-y2)^2 );
%         if dist < space
%             connectsMat(ct,1) = i-1; %For Cpp notation (and .vtk counting)
%             connectsMat(ct,2) = i;   %For Cpp notation (and .vtk counting)
%         ct=ct+1;
%         end
%     elseif i==N
%         x1=xLag(N); x2=xLag(1);
%         y1=yLag(N); y2=yLag(1);
%         dist = sqrt( (x1-x2)^2 + (y1-y2)^2 );
%         if dist < space
%             connectsMat(ct,1) = N-1; %For Cpp notation (and .vtk counting)
%             connectsMat(ct,2) = 0;   %For Cpp notation (and .vtk counting)
%         ct=ct+1;
%         end
%     end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives appropriate string number for filename in printing the
% .vtk files.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function strNUM = give_String_Number_For_VTK(num)

%num: # of file to be printed

if num < 10
    strNUM = ['000' num2str(num)];
elseif num < 100
    strNUM = ['00' num2str(num)];
elseif num<1000
    strNUM = ['0' num2str(num)];
else
    strNUM = num2str(num);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes vorticity from two matrices, U and V, where each
% matrix is the discretized field of velocity values either for x or y,
% respectively.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vort = give_Me_Vorticity(U,V,dx,dy)

% w = ( dv/dx - du/dy )\hat{z}

%Compute dv/dx using central differencing! (maintains periodicity)
dvdx = D(V,dx,'x');

%Compute du/dy using central differencing! (maintains periodicity)
dudy = D(U,dy,'y');

%Compute vorticity
vort = ( dvdx - dudy );

%Take transpose so all lines up
vort = vort';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes vorticity from two matrices, U and V, where each
% matrix is the discretized field of velocity values either for x or y,
% respectively.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uMag = give_Me_Magnitude_Velocity(U,V)

% U: x-directed velocity
% V: y-directed velocity

% Compute magnitude of velocity
uMag = ( U.^2 + V.^2 ).^(1/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints matrix vector data to vtk formated file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function savevtk_points_connects( X, filename, vectorName,connectsMat)

%X is matrix of size Nx3

N = length( X(:,1) );
Nc = length( connectsMat(:,1) );

%TRY PRINTING THEM AS UNSTRUCTURED_GRID
file = fopen (filename, 'w');
fprintf(file, '# vtk DataFile Version 2.0\n');
fprintf(file, [vectorName '\n']);
fprintf(file, 'ASCII\n');
fprintf(file, 'DATASET UNSTRUCTURED_GRID\n\n');
%
fprintf(file, 'POINTS %i float\n', N);
for i=1:N
    fprintf(file, '%.15e %.15e %.15e\n', X(i,1),X(i,2),X(i,3));
end
fprintf(file,'\n');
%
fprintf(file,'CELLS %i %i\n',Nc,3*Nc); %First: # of "Cells", Second: Total # of info inputed following
for s=1:Nc
    fprintf(file,'%i %i %i\n',2, connectsMat(s,1), connectsMat(s,2) );
end
fprintf(file,'\n');
%
fprintf(file,'CELL_TYPES %i\n',Nc); % N = # of "Cells"
for i=1:Nc
   fprintf(file,'3 ');
end
fprintf(file,'\n');
fclose(file);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints matrix vector data to vtk formated file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function savevtk_points( X, filename, vectorName)

%X is matrix of size Nx3

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

%TRY PRINTING THEM AS UNSTRUCTURED_GRID
file = fopen (filename, 'w');
fprintf(file, '# vtk DataFile Version 2.0\n');
fprintf(file, [vectorName '\n']);
fprintf(file, 'ASCII\n');
fprintf(file, 'DATASET UNSTRUCTURED_GRID\n\n');
%
fprintf(file, 'POINTS %i float\n', N);
for i=1:N
    fprintf(file, '%.15e %.15e %.15e\n', X(i,1),X(i,2),X(i,3));
end
fprintf(file,'\n');
%
fprintf(file,'CELLS %i %i\n',N,2*N); %First: # of "Cells", Second: Total # of info inputed following
for s=0:N-1
    fprintf(file,'%i %i\n',1,s);
end
fprintf(file,'\n');
%
fprintf(file,'CELL_TYPES %i\n',N); % N = # of "Cells"
for i=1:N
   fprintf(file,'1 ');
end
fprintf(file,'\n');
fclose(file);

%TRY PRINTING THEM AS POLYGONAL DATA
% file = fopen (filename, 'w');
% fprintf(file, '# vtk DataFile Version 2.0\n');
% fprintf(file, [vectorName '\n']);
% fprintf(file, 'ASCII\n');
% fprintf(file, 'DATASET STRUCTURED_GRID\n');
% fprintf(file, 'DIMENSIONS 64 1 1\n');
% fprintf(file, 'POINTS %i float\n', N);
% for i=1:N
% fprintf(file, '%.15e %.15e %.15e\n', X(i,1),X(i,2),X(i,3));
% end
% fprintf(file,'1.1 1.1 0\n');
% fprintf(file,'CELL_DATA 1\n');
% fprintf(file,'POINT_DATA %u \n',N);
% fprintf(file,'FIELD FieldData 1\n');
% fprintf(file,'nodal 1 %i float\n',N);
% fprintf(file,'0 1 1.1 2\n');
% fprintf(file,'SCALARS nodal float\n');
% fprintf(file,['SCALARS ' vectorName ' float 1 \n']);
% fprintf(file,'LOOKUP_TABLE default\n');

% TRY PRINTING THEM AS POINTS
% file = fopen (filename, 'w');
% fprintf(file, '# vtk DataFile Version 2.0\n');
% fprintf(file, 'Cube example\n');
% fprintf(file, 'ASCII\n');
% fprintf(file, 'DATASET UNSTRUCTURED_GRID\n');
% fprintf(file, 'POINTS %i float\n', N);
% for i=1:N
% fprintf(file, '%.15e %.15e %.15e\n', X(i,1),X(i,2),X(i,3));
% end
% fprintf(file,'POINT_DATA %u \n',N);
% fprintf(file,['SCALARS ' vectorName ' float 1 \n']);
% fprintf(file,'LOOKUP_TABLE default\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints matrix vector data to vtk formated file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function savevtk_vector(X, Y, filename, vectorName,dx,dy)
%  savevtkvector Save a 3-D vector array in VTK format
%  savevtkvector(X,Y,Z,filename) saves a 3-D vector of any size to
%  filename in VTK format. X, Y and Z should be arrays of the same
%  size, each storing speeds in the a single Cartesian directions.
    if (size(X) ~= size(Y))
        fprint('Error: velocity arrays of unequal size\n'); return;
    end
    [nx, ny, nz] = size(X);
    fid = fopen(filename, 'wt');
    fprintf(fid, '# vtk DataFile Version 2.0\n');
    fprintf(fid, 'Comment goes here\n');
    fprintf(fid, 'ASCII\n');
    fprintf(fid, '\n');
    fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
    fprintf(fid, 'DIMENSIONS    %d   %d   %d\n', nx, ny, nz);
    fprintf(fid, '\n');
    fprintf(fid, 'ORIGIN    0.000   0.000   0.000\n');
    %fprintf(fid, 'SPACING   1.000   1.000   1.000\n'); if want [1,32]x[1,32] rather than [0,Lx]x[0,Ly]
    fprintf(fid, ['SPACING   ' num2str(dx) ' '  num2str(dy) '   1.000\n']);
    fprintf(fid, '\n');
    fprintf(fid, 'POINT_DATA   %d\n', nx*ny);
    fprintf(fid, ['VECTORS ' vectorName ' double\n']);
    fprintf(fid, '\n');
    for a=1:nz
        for b=1:ny
            for c=1:nx
                fprintf(fid, '%f ', X(c,b,1));
                fprintf(fid, '%f ', Y(c,b,1));
                fprintf(fid, '%f ', 0);
            end
            fprintf(fid, '\n');
        end
    end
    fclose(fid);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints scalar matrix to vtk formated file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function savevtk_scalar(array, filename, colorMap,dx,dy)
%  savevtk Save a 3-D scalar array in VTK format.
%  savevtk(array, filename) saves a 3-D array of any size to
%  filename in VTK format.
    [nx, ny, nz] = size(array);
    fid = fopen(filename, 'wt');
    fprintf(fid, '# vtk DataFile Version 2.0\n');
    fprintf(fid, 'Comment goes here\n');
    fprintf(fid, 'ASCII\n');
    fprintf(fid, '\n');
    fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
    fprintf(fid, 'DIMENSIONS    %d   %d   %d\n', nx, ny, nz);
    fprintf(fid, '\n');
    fprintf(fid, 'ORIGIN    0.000   0.000   0.000\n');
    %fprintf(fid, 'SPACING   1.000   1.000   1.000\n'); if want [1,32]x[1,32] rather than [0,Lx]x[0,Ly]
    fprintf(fid, ['SPACING   ' num2str(dx) ' '   num2str(dy) '   1.000\n']);
    fprintf(fid, '\n');
    fprintf(fid, 'POINT_DATA   %d\n', nx*ny*nz);
    fprintf(fid, ['SCALARS ' colorMap ' double\n']);
    fprintf(fid, 'LOOKUP_TABLE default\n');
    fprintf(fid, '\n');
    for a=1:nz
        for b=1:ny
            for c=1:nx
                fprintf(fid, '%d ', array(c,b,a));
            end
            fprintf(fid, '\n');
        end
    end
    fclose(fid);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Lagrangian pt data w/ associated scalar to vtk formated file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function savevtk_points_with_scalar_data( X, scalarArray, filename, vectorName)

%X is matrix of size Nx3

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

%TRY PRINTING THEM AS UNSTRUCTURED_GRID
file = fopen (filename, 'w');
fprintf(file, '# vtk DataFile Version 2.0\n');
fprintf(file, [vectorName '\n']);
fprintf(file, 'ASCII\n');
fprintf(file, 'DATASET UNSTRUCTURED_GRID\n\n');
%
fprintf(file, 'POINTS %i float\n', N);
for i=1:N
    fprintf(file, '%.15e %.15e %.15e\n', X(i,1),X(i,2),X(i,3));
end
fprintf(file,'\n');
%
fprintf(file, 'POINT_DATA   %d\n', N);
fprintf(file, ['SCALARS ' vectorName ' double\n']);
fprintf(file, 'LOOKUP_TABLE default\n');
fprintf(file, '\n');
    for i=1:N
        fprintf(file, '%d ', scalarArray(i,1));
        fprintf(file, '\n');
    end

fclose(file);

%     for a=1:nz
%         for b=1:ny
%             for c=1:nx
%                 fprintf(fid, '%d ', array(c,b,a));
%             end
%             fprintf(fid, '\n');
%         end
%     end

% fprintf(file,'CELLS %i %i\n',N,2*N); %First: # of "Cells", Second: Total # of info inputed following
% for s=0:N-1
%     fprintf(file,'%i %i\n',1,s);
% end
% fprintf(file,'\n');
% %
% fprintf(file,'CELL_TYPES %i\n',N); % N = # of "Cells"
% for i=1:N
%    fprintf(file,'1 ');
% end
% fprintf(file,'\n');
% fclose(file);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of vertex pts and all the vertex pts from the
%           .vertex file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [N,xLag,yLag] = read_Vertex_Points(struct_name)

filename = [struct_name '.vertex'];  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

vertices = C{1};    %Stores all read in data in vertices (N+1,2) array

N = vertices(1,1);  % # of Lagrangian Pts
xLag = zeros(N,1);  % Initialize storage for Lagrangian Pts.
yLag = xLag;        % Initialize storage for Lagrangian Pts.

for i=1:N
   xLag(i,1) = vertices(i+1,1); %Stores x-values of Lagrangian Mesh
   yLag(i,1) = vertices(i+1,2); %Stores y-values of Lagrangian Mesh

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of tracer pts and all the tracer pts from the
%           .tracer file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [N,xLag,yLag] = read_Tracer_Points(struct_name)

filename = [struct_name '.tracer'];  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

tracers = C{1};    %Stores all read in data in vertices (N+1,2) array

N = tracers(1,1);  % # of Lagrangian Pts
xLag = zeros(N,1);  % Initialize storage for Lagrangian Pts.
yLag = xLag;        % Initialize storage for Lagrangian Pts.

for i=1:N
   xLag(i,1) = tracers(i+1,1); %Stores x-values of Tracer Lagrangian Mesh
   yLag(i,1) = tracers(i+1,2); %Stores y-values of Tracer Lagrangian Mesh

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the diffusion coefficient and initial concentration, C
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,kDiff] = read_In_Concentration_Info(struct_name,N)

filename = [struct_name '.concentration'];  %Name of file to read in

strstr = '%f';
for i=1:N-1
   strstr = [strstr ' %f'];
end

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,strstr,'CollectOutput',1);

fclose(fileID);        %Close the data file.

con_info = C{1};    %Stores all read in data

%Store all elements on .concentration file
kDiff = con_info(1,1);     %coefficient of diffusion
C = con_info(2:end,1:end); %initial concentration

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of springs and all MASTER NODEs, SLAVE NODEs,
%           spring STIFFNESSES, spring RESTING LENGTHS, spring linearities
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs = read_Spring_Points(struct_name)

filename = [struct_name '.spring'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

spring_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .spring file into a matrix starting w/ 2nd row of read in data to 4th col.
springs = spring_info(2:end,1:4);

%Stores last column to check for non-linearities
NLspr = spring_info(2:end,5);

%Finds indices of non-linearities
inds = isnan(NLspr);
springs2 = zeros(length(inds),1);

%Stores linearity coefficients
for i=1:length(inds)
    if inds(i) == 1
        springs2(i) = 1;
    else
        springs2(i) = NLspr(i);
    end
end

springs = [springs springs2];

clear springs2 inds NPspr spring_info;

%springs: col 1: starting spring pt (by lag. discretization)
%         col 2: ending spring pt. (by lag. discretization)
%         col 3: spring stiffness
%         col 4: spring resting lengths
%         col 5: linearity of spring (1=linear, > 1 non-linear spring)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of TARGET PTS, TARGET-PT-NODEs, and their
%           Target-STIFFNESSES
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function masses = read_Mass_Points(struct_name)

filename = [struct_name '.mass'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f','CollectOutput',1);

fclose(fileID);      %Close the data file.

mass_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .mass file into a matrix starting w/ 2nd row of read in data.
masses = mass_info(2:end,1:3);

%masses:  col 1: Lag Pt. ID w/ Associated Mass Pt.
%         col 2: "Mass-Spring" stiffness Parameter
%         col 3: Mass Value Parameter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of TARGET PTS, TARGET-PT-NODEs, and their
%           Target-STIFFNESSES
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = read_Target_Points(struct_name)

filename = [struct_name '.target'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

targets_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .spring file into a matrix starting w/ 2nd row of read in data.
targets = targets_info(2:end,1:2);

%targets: col 1: Lag Pt. ID w/ Associated Target Pt.
%         col 2: target STIFFNESSES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of POROUS PTS, POROUS-PT-NODEs, and their
%           POROUSITY-COEFFICIENTS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function porosity = read_Porous_Points(struct_name)

filename = [struct_name '.porous'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

porous_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .spring file into a matrix starting w/ 2nd row of read in data.
porosity = porous_info(2:end,1:3);

%porous:  col 1: Lag Pt. ID w/ Associated Porous Pt.
%         col 2: Porosity coefficient

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of beams and all 1st Pt, MIDDLE Pt, and 3rd Pt
%           beam STIFFNESSES, and CURVATURE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function beams = read_Beam_Points(struct_name)

filename = [struct_name '.beam'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f','CollectOutput',1);

fclose(fileID);      %Close the data file.

beam_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .beam file into a matrix starting w/ 2nd row of read in data.
beams = beam_info(2:end,1:5);

    %beams:      col 1: 1ST PT.
    %            col 2: MIDDLE PT. (where force is exerted)
    %            col 3: 3RD PT.
    %            col 4: beam stiffness
    %            col 5: curavture

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of NON-INVARIANT beams and all 1st Pt, MIDDLE Pt,
%           and 3rd Pt beam STIFFNESSES, and CURVATURE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function beams = read_nonInv_Beam_Points(struct_name)

filename = [struct_name '.nonInv_beam'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f %f','CollectOutput',1);

fclose(fileID);      %Close the data file.

beam_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .beam file into a matrix starting w/ 2nd row of read in data.
beams = beam_info(2:end,1:6);

    %beams:      col 1: 1ST PT.
    %            col 2: MIDDLE PT. (where force is exerted)
    %            col 3: 3RD PT.
    %            col 4: beam stiffness
    %            col 5: curavture

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of beams and all 1st Pt, MIDDLE Pt, and 3rd Pt
%           beam STIFFNESSES, and CURVATURE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dSprings = read_Damped_Spring_Points(struct_name)

filename = [struct_name '.d_spring'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f','CollectOutput',1);

fclose(fileID);      %Close the data file.

dSprings_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .d_springs file into a matrix starting w/ 2nd row of read in data.
dSprings = dSprings_info(2:end,1:5);

% dSprings: col 1: starting spring pt (by lag. discretization)
%           col 2: ending spring pt. (by lag. discretization)
%           col 3: spring stiffness
%           col 4: spring resting lengths
%           col 5: damping coefficient

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of muscles and all MASTER NODEs, SLAVE NODEs,
%           length for max. muscle tension, muscle constant, hill
%           parameters (a and b), and Force-Max
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function muscles = read_Muscle_Points(struct_name)

filename = [struct_name '.muscle'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

muscle_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .muscle file into a matrix starting w/ 2nd row of read in data.
muscles = muscle_info(2:end,1:7);

%muscles: col 1: MASTER NODE (by lag. discretization)
%         col 2: SLAVE NODE (by lag. discretization)
%         col 3: length for max. muscle tension
%         col 4: muscle constant
%         col 5: hill parameter, a
%         col 6: hill parameters, b
%         col 7: force maximum!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of muscles and all MASTER NODEs, SLAVE NODEs,
%           length for max. muscle tension, muscle constant, hill
%           parameters (a and b), and Force-Max
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function muscles = read_Hill_3Muscle_Points(struct_name)

filename = [struct_name '.3_muscle'];  %Name of file to read in

fileID = fopen(filename);

    % Read in the file, use 'CollectOutput' to gather all similar data together
    % and 'CommentStyle' to to end and be able to skip lines in file.
    C = textscan(fileID,'%f %f %f %f %f %f %f %f %f','CollectOutput',1);

fclose(fileID);        %Close the data file.

muscle_info = C{1};    %Stores all read in data in vertices (N+1,2) array

%Store all elements on .muscle file into a matrix starting w/ 2nd row of read in data.
muscles = muscle_info(2:end,1:9);

%muscles: col 1: MASTER NODE (by lag. discretization)
%         col 2: SLAVE NODE (by lag. discretization)
%         col 3: length for max. muscle tension
%         col 4: muscle constant
%         col 5: hill parameter, a
%         col 6: hill parameters, b
%         col 7: force maximum!
%         col 8: NL Spring stiffness, kSpr
%         col 9: NL Spring deg. of non-linearity, alpha

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: READS IN ALL THE DATA FOR THE USER-DEFINED FORCE FUNCTION!!!
%           NOTE: DOES NOT SPECIFY HOW MANY PARAMETERS THERE ARE.
%           NOTE: COMPLETELY USER-DEFINED
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function force_general = read_General_Forcing_Function(struct_name)

filename = [struct_name '.user_force'];  %Name of file to read in

% Imports all the data into a data structure
gen_force_info = importdata(filename,' ',1);

% Save data to new variable, force_general
force_general = gen_force_info.data;

author: nickabattista


text:
function edge = clipLineRect(line, rect)
%CLIPLINERECT clip a line with a polygon.
%
%   EDGE = clipLineRect(LINE, RECT);
%   LINE: line in the form [x0 y0 dx dy]
%   RECT: a rectangle in the form [xr yr wr hr] (xr and yr: coordinate of
%   first point, wr and hr are width and height of rectangle)
%
%   Deprecated: use function clipLine instead

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2006-07-24
% Copyright 2006 INRA - Cepia Software Platform

if size(line, 1)==1
    line = repmat(line, size(rect, 1), 1);
elseif size(rect, 1)==1
    rect = repmat(rect, size(line, 1), 1);
elseif size(line, 1) ~= size(rect, 1)
    error('bad sizes for input');
end

edge = zeros(size(line, 1), 4);
for i=1:size(line, 1)
    x = rect(i, 1); y = rect(i, 2); w = rect(i, 3); h = rect(i, 4);

	% intersection with axis : x=xmin
	px1 = intersectLineEdge(line(i,:), [x y x+w y]);
	px2 = intersectLineEdge(line(i,:), [x+w y x+w y+h]);
	py1 = intersectLineEdge(line(i,:), [x+w y+h x y+h]);
	py2 = intersectLineEdge(line(i,:), [x y+h x y]);

	% sort points along the x coordinate, and  draw a line between
	% the two in the middle
	points = sortrows([px1 ; px2 ; py1 ; py2], 1);
	if points(2,1)>=x-1e-14 && points(2,1)<=x+w+1e-14
        if isfinite(points(3,1))
            edge(i, 1:4) = [points(2,:) points(3,:)];
        else
            edge(i, 1:4) = [points(1,:) points(2,:)];
        end
    else
        % line outside the rectangle
        edge(i, 1:4) = [0 0 0 0];
	end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function mps = mph2mps(mph)
%MPH2MPS Convert speed from miles per hour to meters per second
%
%  mps = MPH2MPS(mph) convert speeds from miles per hour to meters per
%   second.
%
%  See also MPH2KMPH, MPH2FTPS, MPH2KTS, MPS2MPH.

% Jonathan Sullivan
% Original: May 2011
% jonathan.sullivan@ll.mit.edu

mps = mph*0.44704;