%% Design of magnetic spring for vertical vibration harvester based on two cuboid magnets (parallel computing)
% Author: Isaac Royo Silvestre (UPNA - Public University of Navarra, Spain)
% Contact information: isaac.royo@unavarra.es
% Version: 2/09/2023

% This is the parallel computing version of harvester_op.m
% Given the desired characteristics (fres and acceleration), a search range 
% (min max dimensions), and a materials list, this script produces a list 
% of magnetic spring models meeting the desired conditions,
% performing an optimization.
% The search/sweep is done in discrete steps, for magnets similar to commercial magnets.
% Magnets min max dimensions should be chosen based on commercial magnets limits,
% practical manufacturing issues, or in order to shrink the search range decreasing
% the run time.
% Magnetic springs are comprised of one fixed magnet and one moving magnet, both with
% cuboid shapes lengthxlengthxheight. It is a vertical system for kinetic harvesters.
% Optimization consists of: For every possible combination of moving magnet and fixed
% magnet lenght and material, minimum fixed magnet height that provies fres is found.
% Thus, the solution list may contain thousands of different solutions,
% deciding which is the best solution is a task for the user, based on their own criteria,
% by means outside the scope of this script.
% Script for MATLAB 2017 or later
% Notes: Implemented magnetic Force equations by G.Akoun and J.P.Yonnet (1984)
%       This is the parallel computing version of harvester_op.m (v
%       6/05/2023), managed by the Matlab Parallel Computing Toolbox.
%       The implementation used the parfor matlab command (the most 
%       transparent but less powerful way to do parallel processing in
%       matlab). The first time a parallel task is run, it may require
%       a few initial seconds while the parallel workers pool is
%       initialized.
%
% Inputs:
% Many constants to be modified by the user in this script (S.I. units)
%
% Outputs:
% sols: a matrix of solutions, 9 columns, undetermined number of rows.
%		Each row contains a magnetic spring model or solution.
%		Column order: fixed magnet legnth, fixed magnet height, moving magnet
%		lenght, moving magnet height, fixed magnet material index, moving magnet material index,
%		z0 (magnets center-to-center equilibrium position distance), kmag (stiffness), fres (Hz)
%		SI units.
%       Notes: This matrix may contain additional columns with "performance
%       indicators" calculated postprocess. This variable can be accesed in the
%		matlab Workspace.
% sols.mat: Saved as binary matlab file in the execution folder.
%       sols.mat is a struct that contains multiple variables, the solution matrix
%       sols, and also data required to interpret
%       the solution or perform additional calculations (the material lists
%       J, magnetro, as well as dampingr, Zp, A0)
% sols_pi.mat: As sols.mat, but the sols matrix includes some performance
%       indicators already calculated.
% harvester_optimizer_log.txt: A text log


%% INPUTS
% Units in SI (m,T,kg)
% Notes: The material index in the index in the J and magnetro vectors
Hmin=0.004; %minimum height of the moving magnet (in practice it wont move if H too small)
Lmin=0.004; %minimum length of the moving magnet
Vmax=0.000003000; %max volume of the moving magnet
HLratiomax=4.5; %max height / length ratio of the moving magnet (this determines the max dimensions of the moving magnet)
HLratiomin=1.2; %min height / length ratio of the moving magnet (limited in practice by friction and magnet rotation possibility)
J=[0.4,1.3]; %possible magnet magnetizations (T), an array with 1 or more values. Use positive numbers.
magnetro=[4900,7500]; %possible magnet material densities (kg/m^3), an array coherent with the J vector (material 1 is index 1 at ro and J). Use positive numbers.
lmin=0.100; %minimum length of the fixed magnet
lmax=0.300; %maximum length of the fixed magnet
hmin=0.001; %minimum height of the fixed magnet
hmax=0.030; %maximum height of the fixed magnet
fres=1; %desired resonant frequency (Hz)
errorf=0.1; %absolute error for fres search (Hz)
dampingr=0.1; %stimated damping ratio of the inertial mass (adimensional, >0)
Zp=0.001; %dead practical space, distance between surfaces of the magnets even at peak displacement (e.g. 1 mm to account the plastic of the frame). Positive value.
A0=0.04; %acceleration amplitude 0-peak (m/s^2) that should be supported, or more. (this increases the dead space so, for A0, magnets cant collide). Positive value.
dL=0.001; %moving magnet lenght search step
dH=0.001; %moving magnet height search step
dh=0.001; %fixed magnet height search step
dl=0.01; %fixed magnet lenght search step
zerrorallowed=0.0001; %absolute error for z0 equilibrium position search error


%% CHECK INPUTS

%Check J ro vectors coherence
if( length(magnetro) ~= length(J) )
    'Error: Materials badly defined, sizes of J and magnetro should be equal'
    return %abort
end

%Correct user inputs
J=abs(J);
magnetro=abs(magnetro);
Zp=abs(Zp);
A0=abs(A0);
dL=abs(dL);
dH=abs(dH);
dh=abs(dh);
dl=abs(dl);
zerrorallowed=abs(zerrorallowed);
dampingr=abs(dampingr);
if(dampingr<=0)
    'Error: Damping must be >0'
    return; %abort
end
errorf=abs(errorf);
fres=abs(fres);
Hmin=abs(Hmin);
Lmin=abs(Lmin);
Vmax=abs(Vmax);
HLratiomax=abs(HLratiomax);
HLratiomin=abs(HLratiomin);
lmin=abs(lmin);
lmax=abs(lmax);
hmin=abs(hmin);
hmax=abs(hmax);
if(lmin>lmax || hmin>hmax || HLratiomin>HLratiomax)
    'Warning: A dimension minimum exceeds its own max'
end


%% ESTIMATE SOME PARAMETERS AND REQUIRED ITERATIONS

%calculated variables
Lmax=min([sqrt(Vmax/Hmin), nthroot(Vmax/HLratiomin,3)]); %inertial mass max side lenght
Y0=A0/(2*pi*fres)^2; %max displacement amplitude, happens for the allowed acceleration at fres
Hmax=min([Vmax/Lmin/Lmin; nthroot(Vmax*(HLratiomax^2),3) ]); %max height for a volume is given for the smaller lenght, but it is limited by the H/L ratio

%Approximated number of required iterations
itmaxAprox=(Lmax-Lmin)/dL*(Hmax-Hmin)/dH*length(J)*length(J)*(lmax-lmin)/dl;

%Exact calculation of the number of required iterations
itmax=0;
Ltemp=Lmin;
while Ltemp<=Lmax+10*eps(Lmax) %eps floating point error
    Hmax=min([Vmax/Ltemp/Ltemp HLratiomax*Ltemp]);
    H=Hmin;
    while H<=Hmax+10*eps(Hmax) 
        %loop for different inertial mass dimensions, at this point
        %a particular inertial mass with a particular H,L is defined
        
        %first, check that the magnet has the allowed H/L ratio
        hlratio=H/Ltemp;
        if(hlratio<=HLratiomax+10*eps(HLratiomax) && hlratio>=HLratiomin-10*eps(HLratiomin))
            itmax=itmax+1;
            
        end %end if h/l ratio check
     
        H=H+dH;
    end
     
    Ltemp=Ltemp+dL;
end
itmax=ceil(itmax*length(J)*length(J)*(lmax-lmin)/dl);

%screen messages
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HARVESTER_OPTIMIZER:\n' );
  fprintf ( 1, '  A script to design magnetic springs for vertical vibration harvesters based on two cuboid magnets.\n' );
  fprintf ( 1, '  (parallel computing version).\n' );
  fprintf ( 1, '  This code was written by Isaac Royo Silvestre (UPNA)\n' );
  fprintf ( 1, '  \n' );
  fprintf ( 1, 'Number of main iterations (find h for a given combination of other parameters) will be less than:\n' );
  fprintf ( 1,'%g',itmaxAprox);
  fprintf ( 1, '  \n' );
  fprintf ( 1, 'and exactly (note that not every iteration provides a particular solution):\n' );
  fprintf ( 1,'%g',itmax);
  fprintf ( 1, '  \n' );



%% MAIN PROCESS
tic %so run time can be registered
L=Lmin;

%First, in order to use a parallel loop "parfor", we save the parameters of
%each main interation (in the sequential script there were multiple nested 
%while loops, something that isn't inmediately translated into a parfor)
itparams=[];
while L<=Lmax+10*eps(Lmax)
    Hmax=min([Vmax/L/L HLratiomax*L]);
    H=Hmin;
    while H<=Hmax+10*eps(Hmax) 
        %loop for different inertial mass dimensions, at this point
        %a particular inertial mass with a particular H,L is defined
        
        %first, check that the magnet has the allowed H/L ratio
        hlratio=H/L;
        if(hlratio<=HLratiomax+10*eps(HLratiomax) && hlratio>=HLratiomin-10*eps(HLratiomin))
                       
            %The inertial mass (moving magnet) is made of one of many possible materials
            materialmm=1; %material index
            while materialmm<=length(J)

                %The fixed magnet is also made of one of many possible materials
                materialfm=1; %material index. Force depends of Jim*Jfm, but spring fres is modified by the inertial mass weight (material density), thus, there are no redundant material combinations
                while materialfm<=length(J)
                    % This will be a main iteration to parfor  (moving magnet decided, fixedmagnet material decided)  
                    itparams=[itparams;
                            [L, H, materialmm, materialfm] ]; % the whole itparams matrix is COPIED to each parallel process, so we keep it as compact as possible
                    materialfm=materialfm+1;
                end
                       
                materialmm=materialmm+1;
            end         
        else %else of the h/l ratio check if
			%do nothing
        end %end of the h/l ratio check if 
        H=H+dH;
    end   
    L=L+dL;
end

%fprintf ( 1, '  \n Progress (iteration/max_iterations): \n' ); %this progress feedback not possible with the current implementation

%Second, now we can use a parfor loop to parallelize the processes
sols=[]; %matrix that holds valid solutions
it=0; %number of main iterations processed, not used because parfor restrictions don't allow updating the results matrix and the iteration counter at the same time
parfor parind=1:1:size(itparams,1)
    L=itparams(parind,1);
    H=itparams(parind,2);
    materialmm=itparams(parind,3);
    materialfm=itparams(parind,4);
    
    %The inertial mass is defined, compute its mass (kg)
    imass=magnetro(materialmm)*L*L*H;
    %Calculate the desired kmag
    kmag0=imass*(2*pi*fres)^2; %k of the virtual spring depends on fres mass etc..
    %Compute the absolute error for kmag, knowing the allowed absolute error for fress (errorf)
    errork=imass*(2*pi*(fres+errorf))^2-kmag0; %in some sources imass*4*pi^2*(errorf^2-2*errorf*fres)
    
    ZA=Y0/2/dampingr; %0-peak max movement amplitude at resonance, according to harmonic movement equation, simplified for f=fres

    subit=0; %to count the number of iterations as would be counted in the sequential script (not the parallel iterations)

    %Now, inertial mass is already defined, as well as
    %the fixed magnet magnetization...
    %We have to check different fixed magnet sizes
    l=lmin; %begin with minimum lenght of the fixed magnet
    complete=false;
    hguess=0; %initial guess for the fixed magnet height, 0 (no guess) for the first iteration
    while l<=lmax && complete==false                   
        %Then, we know everything except the minimum required
        %height for the fixed magnet so the moving magnet
        %is subject to the Force adequate to leviate and
        %do so with the desired kmag...

        %optimize the height of the fixed magnet for this
        %lenght
        Jfm=J(materialfm);
        Jmm=J(materialmm);
        h=optimize_hb(l,L,H,imass*9.8,Jfm,Jmm,Zp,ZA,zerrorallowed,hmax,kmag0,errork,hmin,dh,hguess);
        if(h(1)>0)
            %if a solution is found, record it,
            %remember that h also contains z0 kmag for
            %that solution, to avoid adittional calculations
            sols=[sols;
                    [l,h(1),L,H,materialfm,materialmm,h(2),h(3),1/2/pi*sqrt(h(3)/imass)] ];
            hguess=h(1);
        elseif(h(1)<0)
            hguess=-1*h(1);
        end

        subit=subit+1;

        if( (h(1)>10*eps(0) && h(1)<=hmin+10*eps(hmin)) || (h(1)<-10*eps(0) && abs(h(1))<=hmin+10*eps(hmin)) )
            %if the last found solution was a fixed magnet
            %with small height, the following iterations
            %(with longer magnets) will provide solutions
            %with shorter magnets o no valid solution
            %at all
            complete=true;
            subit=subit+(lmax-l)/dl; %add the number of skipped iterations
        end

        if uint16((it-500*floor(it/500)))==0
            %print the main iteration progress on screen
            fprintf ( 1, 'progress: 500 iterations completed' );
            fprintf ( 1, '\n' );        
        end

        l=l+dl;
    end %end while
    
    fprintf ( 1, 'A parallel process completed: #' );fprintf ( 1,'%u',parind); fprintf ( 1, ' of ' ); fprintf ( 1,'%u',size(itparams,1));
    fprintf ( 1, '\n' );   
    
end %en parfor

toc %show the run time until now

%save solutions to disk (a mat file), and data required to interpret
%the solution or perform additional calculations
save('sols.mat','sols','J','magnetro','dampingr','Zp','A0');

fprintf ( 1, '\n Optimization of the magnetic spring finished \n' );

%postprocess (compute performance indicators and merit figures)
fprintf ( 1, '\n Performance indicators calculation begins, postprocess \n' );
sols=[sols,zeros(size(sols,1),1)]; %new column to te right of the sols matrix, to save the indicators
for i=1:size(sols,1)
	sols(i,end)=compute_vpi(sols(i,:),J,magnetro,0.1,0.5*9.8,350); %for every particular solution, compute and record the Vpi performance indicator
end

%save sols to disk again, but now includes performance indicators
save('sols_pi.mat','sols','J','magnetro','dampingr','Zp','A0');

%end
fprintf ( 1, '\n' );
fprintf ( 1, ' Normal end of execution.\n' );
toc


%% WRITE LOG, DATA IN TEXT FILE
fprintf ( 1, '\n' );
fprintf ( 1, ' Saving log.txt to disk.\n' );
reportname="harvester_optimizer_log.txt";
%header
fid = fopen(reportname,'wt');
fprintf(fid, 'Harvester Optimizer Log\nSummary of results\nLfixedmagnet(m)\tHfixedmagnet(m)\tLmovingmagnet(m)\tHmovingmagnet(m)\tMaterialfixedmagnet(index)\tMaterialmovingmagnet(index)\tz0(m)\tkmag(N/m)\tfres(Hz)\tVpi(V)\n' );
fclose(fid);
for i=1:size(sols,1)
    %fid = fopen(reportname,'at');
    %fprintf(fid, '\n' );
    %fclose(fid);
    dlmwrite(reportname,sols(i,:),'-append','precision','%16.6E','Delimiter','\t');
end



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS

%% optimize_hb(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,errork,hmin,dh,varargin) returns the minimum (optimum) fixed magnet height h required to obtain a magnetic spring with the desired kmag0 (a guess is optional)
% Given all the characteristics of the inertial mass (moving magnet, weight, etc)
% and of the fixed magnet (side, magnetization,...)
% EXCEPT the height of the fixed magnet h,
% this function finds the (aproximated) smallest height the fixed magnet
% should have to provide a "valid" virtual spring with a kmag close to
% the derired stiffness (we assume kmag<0, or Force spatial derivative
% negative)
% See compute_Z0b for notes about signs and assumptions
% This function calls functions compute_hbinterval compute_hbintervaltwice
% when neccesary.
% Inputs:
%  l the side lenght (l) of the origin magnet (fixed), h (height) is the unknown variable
%  L,H the side lenght (L) and height (H) of the target magnet (moving or levitating)
%  Fw the weight of the moving mass, Newtons in SI (this weight must 
%  add the magnet weight and any glued dummy mass weight)
%  J magnetization of the origin magnet
%  Jprima magnetization of the target magnet
%  Zp is the dead space (clearance) to force between magnets faces
%  (might be 0, then a minimal ammount will be used for Zp)
%  ZA is the moving magnet movement amplitude, then the distance between 
%  both magnets must be at least Zp+ZA (ZA can be 0).
%  error is the allowed absolute error (for z0 distance)
%  hmax is the max allowed value for the fixed magnet height h
%  kmag0 is the desired stifness for the magnetic spring (N/m in SI)
%  hmin is the min allowed value for the fixed magnet height h
%  errork is the allowed absolute error for kmag (if 0, by default 
%  100% error is allowed )
%  dh discrete step in the possible values of the height h
%  varargin, an optional input is hguess (an initial guess for h)
% Output: A vector with three numbers [result,z0,kmag0]
%  result is 0 (no valid solution found) or a valid number h
%  for the fixed magnet height.
%  result is -h if a solution was found, but didnt comply with errork and
%  we had no guess (so this h can be used as next guess)
%  z0, kmag will be 0 unless a valid solution was found
function valor=optimize_hb(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,errork,hmin,dh,varargin)
    hguess=0;
    if(length(varargin)>0)
        hguess=varargin{1};
    end
    kmag0=-1*abs(kmag0);
    h=hmin; %minimum height of the fixed magnet
    if(errork==0)
        errork=kmag0;
    end
    if(hguess<=h || hguess==0)
        %if no valid guess, call the usual processs
        %search starting at min height
        sol=compute_hbinterval(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,h,errork,dh);
        
        while (int8(sol(2))==-1) && sol(1)<hmax %this procedure may be slow, but can find solutions in bistable, tristable systems, etc (a general procedure).
            %if a solution that complies with zp zA was not found, keep searching up to hmax
            sol=compute_hbinterval(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,sol(1)+dh,errork,dh);
        end

    else
        %else, pass the guess, to search around the guess
        sol=compute_hbintervaltwice(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,h,errork,dh,hguess);
    end

    valor=[0,0,0];
    if(int8(sol(2))==1)
        valor=[sol(1),sol(3),sol(4)];
    elseif( (hguess<=h || hguess==0) && int8(sol(2))==-2)
        valor=[-1*sol(1),0,0];
    end
    return;
end

%% compute_hbinterval(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,hmin,errork,dh,varargin) returns the minimum (optimum) fixed magnet height h required to obtain a magnetic spring with the desired kmag0, search starts at hmin
% Given all the characteristics of the inertial mass (moving magnet, weight, etc)
% and of the fixed magnet (side, magnetization,...)
% EXCEPT the height of the fixed magnet h,
% this function finds the (aproximated) smallest height the fixed magnet
% should have to provide a "valid" virtual spring with a kmag close to
% the derired stiffness (we assume kmag<0, or Force spatial derivative
% negative)
% See compute_Z0b for notes about signs and assumptions
% Inputs:
%  l the side lenght (l) of the origin magnet (fixed), h (height) is the unknown variable
%  L,H the side lenght (L) and height (H) of the target magnet (moving or levitating)
%  Fw the weight of the moving mass, Newtons in SI (this weight must 
%  add the magnet weight and any glued dummy mass weight)
%  J magnetization of the origin magnet
%  Jprima magnetization of the target magnet
%  Zp is the dead space (clearance) to force between magnets faces
%  (might be 0, then a minimal ammount will be used for Zp)
%  ZA is the moving magnet movement amplitude, then the distance between 
%  both magnets must be at least Zp+ZA (ZA can be 0).
%  error is the allowed absolute error (for z0 distance)
%  hmax is the max allowed value for the fixed magnet height h
%  kmag0 is the desired stifness for the magnetic spring (N/m in SI)
%  hmin is the min allowed value for the fixed magnet height h
%  errork is the allowed absolute error for kmag
%  dh discrete step in the possible values of the height h
% Output: A vector with two numbers
%  Returns [h,0] if no root was found when searching up to h
%  Returns [hsol,1] if root at hsol was found (the nearest) being
%  a valid solution.
%  Returns [h,-1] if a root was found in the search interval up to h,
%  but it is not a valid solution(because displacement is excesive, but complies with kerror)
%  Returns [h,-2] if a root was found in the search interval up to h,
%  but it is not a valid solution.
function valor=compute_hbinterval(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,hmin,errork,dh,varargin)
    kmag0=-1*abs(kmag0); %only negative k is valid
    h=hmin;
    z02=compute_Z0b(l,h,L,H,Fw,J,Jprima,Zp,ZA,error);
    %z01=z02;
    dFb1=-1*(compute_dFmzdzb(l,h,L,H,z02,J,Jprima))-kmag0; %the function which root is searched, compute_dFmzdzb is multiplied by -1 because by default it returns the attraction curve (inverse signs), thus the actual spring stiffness is -compute_dFmzdzb
    dFb2=dFb1;
    signa=sign(dFb1);
    signb=signa; %first iteration, we only know one point, the starting point for h search. Bolzano theorem is applied to an interval [a,b]. Initially a=b.
    %According to Bolzano (bisection), when the target function sign swaps,
    %the root is found.
    while (h+dh)<=hmax && signa==signb
        h=h+dh; %move forward to find a sign swap (h to the right limit)
        z01=z02; %the previous right limit 02 becomes into the left limit 01
        z02=compute_Z0b(l,h,L,H,Fw,J,Jprima,Zp,ZA,error); %calculate the new right limit
        dFb1=dFb2; %the previous right limit becomes into the left limit
        dFb2=-1*(compute_dFmzdzb(l,h,L,H,z02,J,Jprima))-kmag0; %the target function
        signb=sign(dFb2);
    end
    %the loop ends beause we located a root (sign a different to sign b)
    %or we got to hmax without locating an interval with a root.
    valor=[h,0,0,0]; %if no roots, return 0
    %however, if there are roots...
    if(signb ~= signa)
        %return the value h that provides k closest to kmag, if valid
        if(abs(dFb2)<=abs(dFb1) && z02>=(h/2+H/2+Zp+ZA) && abs(dFb2)<=abs(errork)) 
            valor=[h,1,z02,abs(dFb2+kmag0)]; %h in the right limit is closer
        elseif(abs(dFb2)>abs(dFb1) && z01>=(h/2+H/2+Zp+ZA) && abs(dFb1)<=abs(errork))
            valor=[h-dh,1,z01,abs(dFb1+kmag0)]; %return the value h that provides k closest to kmag
        else
            valor=[h,-1,0,0]; %desired kmag is in this invertal, but does not comply with error or hmin. We return a not found code
            if( abs(dFb2)>abs(errork) && abs(dFb1)>abs(errork) )
                 valor=[h,-2,0,0];
            end
        end
    end
    return;
end

%% compute_hbintervaltwice(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,hmin,errork,dh,hguess) returns the minimum (optimum) fixed magnet height h required to obtain a magnetic spring with the desired kmag0, using a user guess
% Given all the characteristics of the inertial mass (moving magnet, weight, etc)
% and of the fixed magnet (side, magnetization,...)
% EXCEPT the height of the fixed magnet h,
% this function finds the (aproximated) smallest height the fixed magnet
% should have to provide a "valid" virtual spring with a kmag close to
% the derired stiffness (we assume kmag<0, or Force spatial derivative
% negative)
% See compute_Z0b for notes about signs and assumptions
% hguess is an initial guess for h, search is done AROUND hguess,
% as long as guess not zero, and guess>hmin
% The search is limited by hmin hmax
% Inputs:
%  See function compute_hbinterval
%  hguess is an initial guess for h, (it often speeds up the process)
% Output: A vector with two numbers.
%  See function compute_hbinterval
function valor=compute_hbintervaltwice(l,L,H,Fw,J,Jprima,Zp,ZA,error,hmax,kmag0,hmin,errork,dh,hguess)
    kmag0=-1*abs(kmag0); %only negative k is valid
    if(hguess>hmax)
        hguess=hmax;
    end
    if(hguess<hmin || hguess==0)
        hguess=hmin;
    end
    h=hguess;
    %Two search processes will be made, step by step,
    %alternating the left search(decrease h from hguess) and the right search
    % (increase h from hguess).
    %First we move to the left (decrease h)
    %since it is called by a process that aims to minimize h
    %when called.
    %the starting point, both search intervals have zero length
    %initially
    z0left2=compute_Z0b(l,h,L,H,Fw,J,Jprima,Zp,ZA,error);
    z0left1=z0left2;
    z0right2=z0left2;
    z0right1=z0right2;
    dFbleft2=-1*(compute_dFmzdzb(l,h,L,H,z0left2,J,Jprima))-kmag0; %the function which root is searched, compute_dFmzdzb is multiplied by -1 because by default it returns the attraction curve (inverse signs), thus the actual spring stiffness is -compute_dFmzdzb
    dFbleft1=dFbleft2;
    dFbright2=dFbleft2;
    dFbright1=dFbright2;
    signlefta=sign(dFbleft2);
    signleftb=signlefta; 
    signrighta=signlefta;
    signrightb=signrighta;
    hleft=h;
    hright=h;
    direction=false; %direction of the current search (it will vary), false left, true right
    %first iteration, we only know one point, the starting point for h search. Bolzano theorem is applied to an interval [a,b]. Initially a=b.
    %According to Bolzano (bisection), when the target function sign swaps,
    %the root is found.
    while ( (hright+dh)<=hmax || (hleft-dh)>=hmin ) && signrighta==signrightb && signlefta==signleftb %while no root is found and h steps in some direction are possible, continue
        if(direction==true)
            if((hright+dh)<=hmax)
                %search to the right
                hright=hright+dh; %move forward to find a sign swap (h to the right limit)
                z0right1=z0right2; %the previous right limit 02 becomes into the left limit 01
                z0right2=compute_Z0b(l,hright,L,H,Fw,J,Jprima,Zp,ZA,error); %calculate the new right limit
                dFbright1=dFbright2; %the previous right limit becomes into the left limit
                dFbright2=-1*(compute_dFmzdzb(l,hright,L,H,z0right2,J,Jprima))-kmag0; %the target function
                signrightb=sign(dFbright2);
            end
        else
            if (hleft-dh)>=hmin
                %(direction==false)
                %search to the left
                hleft=hleft-dh; %move forward to find a sign swap (h to the left limit)
                z0left2=z0left1; %the previous left limit 01 becomes into the right limit 02
                z0left1=compute_Z0b(l,hleft,L,H,Fw,J,Jprima,Zp,ZA,error); %calculate the new left limit
                dFbleft2=dFbleft1; %the previous left limit becomes into the right limit
                dFbleft1=-1*(compute_dFmzdzb(l,hleft,L,H,z0left1,J,Jprima))-kmag0; %the target function
                signlefta=sign(dFbleft1);
            end
        end
        direction=~direction; %the ~ operator swaps the logic value (true to false, false to true)
    end
    %the loop ends beause we located a root (sign a different to sign b)
    %or we got to hmax without locating an interval with a root.
    valor=[hmax,0,0,0]; %if no roots, return 0
    %however, if there are roots...
    if(signrightb ~= signrighta)
        %return the value h that provides k closest to kmag, if valid
        if(abs(dFbright2)<=abs(dFbright1) && z0right2>=(hright/2+H/2+Zp+ZA) && abs(dFbright2)<=abs(errork)) 
            valor=[hright,1,z0right2,abs(dFbright2+kmag0)]; %h in the right limit is closer
        elseif(abs(dFbright2)>abs(dFbright1) && z0right1>=((hright-dh)/2+H/2+Zp+ZA) && abs(dFbright1)<=abs(errork))
            valor=[hright-dh,1,z0right1,abs(dFbright1+kmag0)]; %return the value h that provides k closest to kmag
        else
            valor=[h,-1,0,0]; %desired kmag is in this invertal, but does not comply with error or hmin. We return a not found code
            if( abs(dFbright2)>abs(errork) && abs(dFbright1)>abs(errork) )
                 valor=[h,-2,0,0];
            end
        end
    elseif(signleftb ~= signlefta)
        %return the value h that provides k closest to kmag, if valid
        if(abs(dFbleft2)<=abs(dFbleft1) && z0left2>=((hleft+dh)/2+H/2+Zp+ZA) && abs(dFbleft2)<=abs(errork)) 
            valor=[hleft+dh,1,z0left2,abs(dFbleft2+kmag0)]; %h in the right limit is closer
        elseif(abs(dFbleft2)>abs(dFbleft1) && z0left1>=((hleft)/2+H/2+Zp+ZA) && abs(dFbleft1)<=abs(errork))
            valor=[hleft,1,z0left1,abs(dFbleft1+kmag0)]; %return the value h that provides k closest to kmag
        else
            valor=[h,-1,0,0]; %desired kmag is in this invertal, but does not comply with error or hmin. We return a not found code
            if( abs(dFbleft2)>abs(errork) && abs(dFbleft1)>abs(errork) )
                 valor=[h,-2,0,0];
            end
        end
    end
    return;
end

%% compute_Z0b(l,h,L,H,Fw,J,Jprima,Zp,ZA,error) returns the equilibrium point z0 between two box shaped magnets in repulsion
% These box magnets are defined lenghtxlenghtxheight (two sides equal),
% centered and aligned in the Z axis.
% To find the Z0 point where both magnets levitate.
% This function assumes that the first magnet (l,h) is the fixed magnet placed below
% the other magnet (L,H) which levitates at certain positive z. Levitation below
% the fixed magnet is not allowed, even if rarely possible.
% It is assumed that J, Jprima are positive (no sign) altough
% magnets have opposite magnetizations (repulsion).
% Calculations use a modified bisection method and calls other
% functions to compute magnetic forces (see compute_Fmagzb and compute_Fmzb)
% Inputs:
%  l,h the side lenght (l) and height (h) of the origin magnet (fixed)
%  L,H the side lenght (L) and height (H) of the target magnet (moving or levitating)
%  Fw the weight of the moving mass, Newtons in SI (this weight must 
%  add the magnet weight and any glued dummy mass weight)
%  J magnetization of the origin magnet
%  Jprima magnetization of the target magnet
%  Zp is the dead space (clearance) to force between magnets faces
%  (might be 0, then a minimal ammount will be used for Zp)
%  ZA is the moving magnet movement amplitude, then the distance between 
%  both magnets must be at least Zp+ZA (ZA can be 0).
%  error is the allowed absolute error (for z0 distance)
% Output: Z0 (double float) is the axial distance between magnets centers
%	while in levitation equilibrium (magnetic Force=Weight)
function valor=compute_Z0b(l,h,L,H,Fw,J,Jprima,Zp,ZA,error)
    %z is the axial distance between magnets centers (called gamma in some
    %formulae). Formulae can fail if the magnets collide, thus a 
    %minimum z is required.
    if(Zp<=0)
        Zp=0.0000001;
    end
    if(ZA<0)
        ZA=0;
    end
    zmin=h/2+H/2+Zp+ZA*0.9;
    z=zmin;
    Fa=-1*(compute_Fmzb(l,h,L,H,z,J,Jprima))-Fw;
    dz=max([l,h,L,H]);
    Zmax=zmin+dz*10;
    while Fa<0 && z<=Zmax %F>0 for it to make sense (stable equilibrium)
        %Otherwise, a point where F>0 has to be found
        z=z+max(h/2,H/2);
        Fa=-1*(compute_Fmzb(l,h,L,H,z,J,Jprima))-Fw;
    end
    
    if(Fa<0)
        %If positive forces are not found in a resonable range,
        %return error. Causes can be the moving magnet being heavier
        %than the levitation force, or maybe the desired movement
        %amplitude exceeds the levitation distance.
        valor=false;
        return;
    end
    
    %Otherwise, continue, save the left limit of the interval
    %And search for the right limit of the interval
    %In the right limit, Force becomes negative
    a=z;
    z=z+dz;
    Fb=-1*(compute_Fmzb(l,h,L,H,z,J,Jprima))-Fw;
    while(Fb>0 && z<=Zmax)
        z=z+dz;
        Fb=-1*(compute_Fmzb(l,h,L,H,z,J,Jprima))-Fw;
    end
    
    if(Fb>0)
        %If the right limit is not found in a resonable range,
        %return error, as we can not guarantee the existence of a root
        %that can be found via Bolzano. It is caused by the magnet weight
        %exceding the levitation force, or rare cases.
        valor=false;
        return;
    end
    
    b=z;
    
    %In case the solution is exactly in one of the interval limits
    if(abs(Fa)<0.000001)
        valor=a;
        return;
    end
    if(abs(Fb)<0.000001)
        valor=b;
        return;
    end
    
    %Finally, apply the bisection method
    n=floor( (log(b-a)-log(error))/log(2) -1 )+1; %number of required iterations according to the bisection method theory
    for i=1:n
        z=(b+a)/2;
        F=-1*(compute_Fmzb(l,h,L,H,z,J,Jprima))-Fw;
        if(F>0)
            a=z;
        else
            b=z;
        end
    end
    Z0=(b+a)/2;
    
    valor=Z0;
    return;
end

%% compute_Fmzb(l,h,L,H,gamma,J,Jprima) returns the axial force between two box shaped magnets with magnetization axis Z and two equal sides
% It just calls the function compute_Fmagzb but simplified:
% Now each magnet has dimensions LxLxH (lenght=wide), so
% these are just defined by side lenght and height.
% Furthermore, there is no need to pass half dimensions as in compute_Fmagzb
% However, both magnets now have to be centered in the axial axis
% (00Z), so only the gamma axial distance is required.
% Note that the obtained Force is the force seens by the target magnet,
% also note the magnetizations signs (see notes in compute_Fmagzb)
% Inputs:
%  l,h the side lenght (l) and height (h) of the origin magnet
%  L,H the side lenght (L) and height (H) of the target magnet
%  gamma is the axial distance between magnets centers
%  J magnetization of the origin magnet
%  Jprima magnetization of the target magnet
% Output: The derivative dFmagz/dz (double float), N/m if inputs in SI units
function valor=compute_Fmzb(l,h,L,H,gamma,J,Jprima)

    valor=compute_Fmagzb(l/2,l/2,h/2,L/2,L/2,H/2,0,0,gamma,J,Jprima);
    return;
end

%% compute_Fmagzb(a,b,c,A,B,C,alpha,beta,gamma,J,Jprima) returns the axial force between two box shaped magnets with magnetization axis Z
% The "origin" magnet, which is recommended to be a fixed magnet located in the bottom, has
% dimensions 2a,2b,2c, is centered in the relative coordinates origin
% 0,0,0, (cartesian coordinates x,y,z) and has magnetization J.
% The "target" magnet, which is recommended to be a moving magnet located
% over the origin magnet, has dimensions 2A,2B,2C, its center placed at the distances  
% alpha,beta,gamma from the originmagnet center in x,y,z axis. Its magnetization is
% Jprima.
% This function returns the axial magntic force (Z axis) seen by
% the target magnet. The force experienced by the origin magnet would the same
% with opposite sign (Newton laws).
% Both magnets have magnetization J axis parallel to Z axis. J is measured in Teslas (SI)
% (this is the naming convention used in the reference paper)
% If J, Jprima are positive, magnets are in attraction (Norths upwards). 
% If signs are opposite, magnets in repulsion (J,Jprima signs indicate the 
% magnetization direction)
% Force will have the correct sign according to the magnets signs
% (if J,Jprima>0 and gamma>0, F<0 if magnets attract. The returned Force can be
% multiplied or give the correct signs to J)
% Use S.I. units (m, T) or at least units coherent with the value of
% mu0 4*pi*10^-7 (vacuum magnetic permeability)
% J is often close to the magnet remanent field Br (i.e. 1.3 T for NdFeB magnets)
% See formulae in D.F. Berdy et al. - "Design and optimization 
% of a magnetically sprung block magnet vibration energy harvester", or
% G.Akoun et al. - "3D ANALYTICAL CALCULATION OF THE FORCES EXERTED BETWEEN TWO CUBOIDAL MAGNETS"
% Warning: In Akoun  paper Fx, Fy, are experimentally corroborated, not Fz.
% This equation produces indeterminations when magnets vertex collide.
% Inputs:
%  a,b,c half sides of the origin magnet (again note that we must pass the half of each side)
%  A,B,C half sides of the target magnet
%  alpha,beta,gamma distances betwen magnets centers
%  J, Jprima magnets magnetizations
% Output: Axial force Fz experienced by the target magnet (double float), N if inputs in SI units
function valor=compute_Fmagzb(a,b,c,A,B,C,alpha,beta,gamma,J,Jprima)
   valor=0;
   %This equation has 256 terms based in 64 combinations of 6
   %index, it can be written as nested loops
   for i = 0 : 1
       for j = 0 : 1
           for k = 0 : 1
               for l = 0 : 1 %index el as seen in the reference paper
                   for p = 0 : 1
                       for q = 0 : 1
                           u=alpha+A*(-1)^j-a*(-1)^i;
                           v=beta+B*(-1)^l-b*(-1)^k;
                           w=gamma+C*(-1)^q-c*(-1)^p;
                           r=sqrt(u^2+v^2+w^2);
                           valor=valor+((-1)^(i+j+k+l+p+q))*((-1)*u*w*log(r-u)-v*w*log(r-v)+u*v*atan(u*v/r/w)-r*w);
                       end
                   end
               end
           end
       end
   end
   valor=valor*J*Jprima/16/pi/pi*10^7; %16 from 4*4
   return;
end

%% compute_dFmzdzb(l,h,L,H,gamma,J,Jprima) returns the derivative of the force Fz to dz between two box shaped magnets with axial magnetization and two equal sides
% It just calls the function compute_dFmagzdzb but simplified:
% Now each magnet has dimensions LxLxH (lenght=wide), so
% these are just defined by side lenght and height.
% Furthermore, there is no need to pass half dimensions as in compute_Fmagzb or dFmagzdzb
% However, both magnets now have to be centered in the axial axis
% (00Z), so only the gamma axial distance is required.
% Note that the obtained value is the derivative of the force seen by the target magnet,
% also note the magnetizations signs (see notes in compute_dFmagzdzb)
% Inputs:
%  l,h the side lenght (l) and height (h) of the origin magnet
%  L,H the side lenght (L) and height (H) of the target magnet
%  gamma is the axial distance between magnets centers
%  J magnetization of the origin magnet
%  Jprima magnetization of the target magnet
% Output: The derivative dFmagz/dz (double float) N/m if inputs in SI
function valor=compute_dFmzdzb(l,h,L,H,gamma,J,Jprima)

    valor=compute_dFmagzdzb(l/2,l/2,h/2,L/2,L/2,H/2,0,0,gamma,J,Jprima);
    return;
end


%% compute_dFmagzdzb(a,b,c,A,B,C,alpha,beta,gamma,J,Jprima) returns the derivative of the force Fz to dz between two box shaped magnets with magnetization axis Z
% Function alike compute_Fmagzb but the sum terms are differentiated
% with respect to z (gamma)
% This provides the stiffness of the function, or spring k at that point
% (as k is the slope of the Force graph)
% Note it is the derivative of the Force experienced by the target magnet
% J signs indicate the magnetization directions, if both
% J,Jprima are positive, then magnets are in attraction
% The equation is obtained by differentiating the magnetic Force equation by
% G.Akoun - 3D ANALYTICAL CALCULATION OF THE FORCES EXERTED BETWEEN TWO CUBOIDAL MAGNETS
% Inputs: See function compute_Fmagzb
% Output: The derivative dFmagz/dz (double float), N/m if inputs in SI units
function valor=compute_dFmagzdzb(a,b,c,A,B,C,alpha,beta,gamma,J,Jprima)
   valor=0;
   %This equation is based in 64 combinations of 6
   %index, it can be written as nested loops
   for i = 0 : 1
       for j = 0 : 1
           for k = 0 : 1
               for l = 0 : 1 %index el as seen in the reference paper
                   for p = 0 : 1
                       for q = 0 : 1
                           u=alpha+A*(-1)^j-a*(-1)^i;
                           v=beta+B*(-1)^l-b*(-1)^k;
                           w=gamma+C*(-1)^q-c*(-1)^p;
                           r=sqrt(u^2+v^2+w^2);
                           valor=valor+((-1)^(i+j+k+l+p+q))*((-1)*(u*log(r-u)+u*w^2/r/(r-u))-(v*log(r-v)+v*w^2/r/(r-v))-((u*v)^2*((w^2)/r+r)/((r*w)^2+(u*v)^2))-(w^2)/r-r);
                       end
                   end
               end
           end
       end
   end
   valor=valor*J*Jprima/16/pi/pi*10^7; %16 from 4*4
   return;
end

%% compute_vpi(sol,J,magnetro,dampingr,Acel,N) returns the Vpi performance indicator
% Vpi is an approximated estimation of the 0-peak voltage the device could generate
% in open circuit.
% Inputs:
%  sol is a row vector that contains one particular solution (a harvester model)
%		Column order: fixed magnet legnth, fixed magnet height, moving magnet
%		lenght, moving magnet height, fixed magnet material index, moving magnet material index,
%		z0 (magnets center-to-center equilibrium position distance), kmag (stiffness), fres (Hz)
%		It may contain additional columns that are not used here.
%  J is a vector that contains the magnetization (T) of each possible material
%  magnetro is a vector that contains the density of each possible material (unused)
%  dampingr estimated damping ratio (adimensional)
%  Acel is the 0-peak acceleration used to calculate Vpi (m/s^2)
%  N is the number of turns in the coil
% Output: Vpi (double float) calculated key performance indicator or merit figure, Volts (SI)
function valor=compute_vpi(sol,J,magnetro,dampingr,Acel,N)
	J=J(sol(6)); %only moving magnet J matters
	L=sol(3);
	H=sol(4);
	keturn=J/pi*( 2*L*sqrt(2) - 2*L^4*(2*L^2+8*H^2) / ((L^2+4*H^2)^2*sqrt(2*L^2+4*H^2)) );
	valor=N*keturn*Acel/4/pi/sol(9)/dampingr;
	return;
end
