%CORELOSS Caclulate power loss in a ferrite core.
%   
%   CORELOSS(time, B, alpha, beta, k, suppressFlag)
%
%   time and B are a piecewise linear series of data points
%   time = vector of successive time values corresponding to flux values, in seconds.
%   B = vector of flux densities corresponding to time vector values, in tesla.
%   alpha, beta, k = Steinmetz parameters.  Must be for use with MKS units
%   suppressFlag: 0 to display input plot, loss  and error information in the
%                   command window.
%                 1 to suppress all output and only return the loss value.
%                   If an error occurs, CORELOSS will return -1.
%             
%   CORELOSS will calculate the power loss in your specified ferrite core.
%   The loss value is in W/m^3. 

% Created by Charles Sullivan, Kapil Venkatachalam, Jens Czogalla
% at Thayer School of Engineering, Dartmouth College.
% Modified by Kip Benson 3/03 (error checking on input data and improved
% command line interface).
% Corrected 11/05 C. Sullivan:  fixed minor loop handling bug ("minortime" was mistyped as minorime")

function y = coreloss(tvs, B, alpha, beta, k, suppressFlag)

%tvs is what was called "time" in help info above.
% check data for errors
	error = 'None';
	
	%makes sure times are successive
	for i = 2 : length(tvs)
        if tvs(i-1) >= tvs(i)
            error = 'Time data points must be successive.';
        end
	end
    
    %make sure vectors are same length
    if length(tvs) ~= length(B)
        error = 'Time and flux vectors must have same length';
    end
    
    if B(1) ~= B( length(B) )
        error = 'Since the PWL input is periodic, the first flux value must equal the last';
    end
% done checking for errors

% calculate core loss if there is no error with the input data
if suppressFlag == 0
    
    if strcmp(error, 'None')
        Pcore = 1000*gsepwl(tvs', B, alpha, beta, k);
        disp(strcat('Core Loss: ', num2str(Pcore), ' (W/m^3)'))
        y = Pcore;
        
        %plot the user's input
        pwlPlot = figure('visible', 'on');
        plot(tvs, B,'b.-');
        xlabel('Time (s)');
        ylabel('Flux Density (T)');
        title('Plot of Input Data');
    else
        disp(strcat('Error: ', error))
        y = -1;
    end
    
elseif suppressFlag == 1                    % suppress output
    if strcmp(error, 'None')
        y = 1000*gsepwl(tvs', B, alpha, beta, k);
    else
        y = -1;     % error occured
    end

end


%====================================================================
%to calculate core loss per unit volume using GSE for PWL signal form
%====================================================================
function p=gsepwl(t,B,alpha,beta,k)

            % a is fraction of Bpp used, 1-a is fraction of original
a = 1.0;    % a=1 for iGSE
            % a=0 for GSE

T = t(length(t))-t(1);            %total time of PWL period

ki = k/((2^(beta+1))*pi^(alpha-1)*(0.2761+1.7061/(alpha+1.354)));

[B,t,Bs,ts] = splitloop(B,t);     %split waveform into major and minor loops

pseg(1) = calcseg(t,B,alpha,beta,ki,a);

dt(1) = t(length(t))-t(1);

for j = 1:length(ts)
    pseg(j+1) = calcseg(ts{j},Bs{j},alpha,beta,ki,a);
    tseg = ts{j};
    dt(j+1) = tseg(length(tseg))-tseg(1);
end

p = sum(pseg.*dt)/T;

%==============
% Loop Splitter
%==============
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Loop Spliter%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%By %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Kapil Venkatachalam%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Modified on June 2nd, 2002%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%This program takes a piecewise linear
%current waveform corresponding to a B-H loop 
%and splits it into smaller waveforms 
%which have the same starting and ending value.
%
%Inputs:
%
%'a': B vector corresponding to the B-H loop.
%'b': Time vector of currents corrssponding to the B-H loop. 
%
%Outputs:
%
%'Majorloop': B values corrssponding to the major B-H loop.
%'Majortime': time vector of currents corrssponding to the major B-H loop.
%'Minorloop': B values corrssponding to the different minor B-H loop.
%'Minorloop': time vector of currents corrssponding to the different minor B-H loop.

function [ majorloop, majortime, minorloop, minortime] = splitloop(a,b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                Reshaping the waveforms and identifying the peak point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Identifies the lowest point in 'a'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

e = 0;                          %Constant defined to identify the position of the lowest value.
sa = 1;                         %Counter to check if lowest point is reached.
for d = 1 : length(a)           %Check for the size of vector.
    if (a(d)== min(a)& sa ==1)  %Checking condition if lowest point has been reached.
        e = d;                  %Identifies the lowest point in the waveform.
        sa = sa +1;             %Counter incremented if lowest point is reached.
    end
end

%v is a vector which stores the shifted values of 'a'.
%t is a vector which stores the shifted values of 'b'.

v = a(e);       % adds the lowest point as the first value.
t = 0;          % adds the coresponding time value the first value.


bdiff = diff(b); % adjusts for the times corresponding to the values.
cumdiff = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loop stores all value of 'a' (from the lowest point till the end point) in 'v'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for q = e+1 : length(a)             %'q' counts the remaing vectors before the lowest point is reached.
    v = [v a(q)];                   %Values of 'a' stored in 'v'.
    cumdiff = cumdiff + bdiff(q-1); %Time adjustment
    t = [t cumdiff];                %Values of 'b' stored in 't'.
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loop stores all value of 'a' (from the starting point till the lowest point) in 'v'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for zx = 2 : e                      % Appends the part before the lowest point.
    v = [v a(zx)];                  %Values of 'a' stored in 'v'.
    cumdiff = cumdiff + bdiff(zx-1);%Time adjustment    
    t = [t cumdiff];                %Values of 'b' stored in 't'.
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         Finds the position of the peak value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

z = 0;                          %Constant defined to identify the position of the peak value.
sa = 1;                         %Counter to check if peak point is reached.
for x= 1 : length(v)            %Check for the size of vector.
    if (v(x)== max(v) & sa == 1)%Checking condition if peak has been reached.
        z = x;                  %Identifies the peak point in the waveform.
        sa = sa +1;             %Counter incremented if peak point is reached.
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             End of reshaping the waveforms and identifying the peak point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Defining variables to keep track of values in the vectors
i=2;
j=1;
k=1;

s=cell(1300); % Defining cells for minorloop.
p=cell(1300); % Defining cells for minortime.

%'m' stores the majorloop values extracted from 'v'.
%'n' stores the corrsponding time values extracted from 't'.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Splits the rising portion of the waveform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m(1) = v(1); % Adds the first element of 'v' to'm'.
n(1) = t(1); % Adds the first element of 't' to'n'.


while (i <= z)          %Checks for all values before the peak point
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                     Check to see if the waveform is rising
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if (v(i) >= v(i-1)) %Compares adjascent values of 'v'.
        
        m = [m v(i)];   %Adds the an element of 'v' to'm'.
        n = [n t(i)];   %Adds the an element of 't' to'n'.
        count=1;        %Counter to keep track of end of the rising part.
        
    else
        %Check for minor loop in the risisng part
        s{j,1} = [s{j,1} v(i-1)];   %Adds the an element of 'v' to's'.
        p{j,1} = [p{j,1} t(i-1)];   %Adds the an element of 't' to'p'.
        
% by jens czogalla
        k = k + 1;
        
        % Repeating the process till the minor loop ends
        while (v(i) < max(m))     
            s{j,1} = [s{j,1} v(i)];     %Adds the an element of 'v' to's'.
            p{j,1} = [p{j,1} t(i)];     %Adds the an element of 't' to'p'.
% by jens czogalla            
            k = k + 1;
            
            i=i+1;                      %Increments the counter keeping track of the elements in 'v'.
        end
        
        % Calculating the slope of the rising edge of the minor loop.
        slope = (v(i-1)-v(i))/(t(i-1)-t(i)); 

        %Makes the last element of the minor loop same as the maximum value of 'm'.
        s{j,1} = [s{j,1} max(m)];
        
        %Computes the value of time of the point which was stored last in 's'.
        stemp = ((max(m) - v(i-1)) / slope) + t(i-1); %Calculating the time value.
        p{j,1} = [p{j,1} stemp];                      %Value is stored in p.
        m = [m max(m)];                               %The last point in 'm' is repeated for continuity.
        n = [n stemp];                                %The time value is also repeated.
        count=count+1;                                %Counter is incremented to indicate end of minor loop.
        j=j+1;                                        %Index of 's' is incremented.
% bz jens czogalla        
        k = 1;
        
    end
   if (count <= 1)                                    %Check condition keeping track of end of rising part.
        i=i+1;                                        %Increment counter keeping track of the elements in 'v'.
   else
   end    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of splitting of the rising portion of the waveform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%'m' now stores part the rising part of the major loop.
%'n' stores the corresponding time values of the major loop.
%'s' stores the minor loops in the rising part of the waveform.
%'p' stores the corresponding time values of the minor loops.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Splits the falling portion of the waveform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while (i <= length(v) ) %Checks for all values after the peak point
    if (v(i) <= v(i-1)) %Compares adjascent values of 'v'.
     
        m = [m v(i)];   %Adds the an element of 'v' to'm'.
        n = [n t(i)];   %Adds the an element of 't' to'n'.
        count=1;        %Counter to keep track of end of the rising part.
        
    else
        
        %Check for minor loop in the falling part.
        temp = v(i-1);                          %Temporary variable to store last value in 'v'.
        s{j,1} = [s{j,1} v(i-1)];               %Adds the an element of 'v' to's'. 
        p{j,1} = [p{j,1} t(i-1)];               %Adds the an element of 't' to'p'. 
% by jens czogalla        
        k = k + 1;
               
        while (i <= length(v) & v(i) > temp)    %Compares adjascent values of 'v' for all the remaining 'v' values.
            s{j,1} = [s{j,1} v(i)];             %Adds the an element of 'v' to's'.
            p{j,1} = [p{j,1} t(i)];             %Adds the an element of 't' to'p'.
% by jens czogalla            
            k = k + 1;
            
            i=i+1;                              %Increments the counter keeping track of the elements in 'v'.
        end   
       
        % Repeating the process till the minor loop ends

%         while (i <= length(v))
% by jens czogalla
        while (i <= length(v) & k > 1)
            
           
            % Calculating the slope of the rising edge of the minor loop.
            slope = (v(i-1)-v(i))/(t(i-1)-t(i));
            
            %Makes the first element of the minor loop same as the last element in 'm'.
            s{j,1} = [s{j,1} temp];
            
            %Computes the value of time of the point which was stored last in 's'.
            qtemp = ((temp - v(i-1)) / slope) + t(i-1); %Calculating the time value.
            p{j,1} = [p{j,1} qtemp];                    %Value is stored in p.
            r=1;                                        %Counter to make the pass through the loop only once.
            while ( v(i) ~= temp & r==1)
                m = [m temp];                           %The last point in 'm' is repeated for continuity.
                n = [n qtemp];                          %The time value is also repeated.
                r = r + 1;                              %Counter incremented to indicate the pass has been made.
            end
            count=count+1;                              %Counter is incremented to indicate end of minor loop.
            j=j+1;                                      %Index of 's' is incremented.
% by jens czogalla            
            k = 1;
            
        end
   end
   if (count <= 1)                                      %Check condition keeping track of end of falling part.
       i=i+1;                                           %Increment counter keeping track of the elements in 'v'.
   else
   end    
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of splitting of the falling portion of the waveform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Removal of repetition of points in 'm' and adjusting the values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x= length (m);          %Variable to store the length of 'm'.
majorloop= [m(1)];      %Stores the first element of 'm' in 'majorloop'                 
majortime= [n(1)];      %Stores the first element of 'n' in 'minorloop'      
g = 0;                  %Initializing variable which adjust time for points following the point of repetition.



% by jens czogalla
for h = 2:x
    

    if m(h-1) ~= m(h)   %Checking for repeated points.
         majortime = [majortime n(h)-g];    %Add time value to 'majortime'.
         majorloop = [majorloop m(h)];      %Add corresponding value of 'm' in 'majorloop'.    
    else
        g = g + n(h) -n(h-1);               %Adjusts the time to compensate for repitions.
    end   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           End of removal of repetition of points in 'm' and adjusting the values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Finding the number of minor loops to be used for checking sub loops later
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    
uo = 1;
pt = 0;                         %Variable to keep track of the number of minor loops.
ss = 1;
while (ss == 1 & uo <= size(s))
        % while (s{uo,1} ~= [])  isempty~(X)
        while ~isempty(s{uo,1})
            uo = uo + 1;
        end
        ss = ss + 1;
        pt = uo-1;
 end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        Recursion
%              Checks if any portion of the  split waveform has subloops.
%              If so repeats the above process to make it a single loop.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

minorloop = {};
minortime = {};
for qq = 1: pt

    sinp = s{qq};    %flux
    pinp = p{qq};    %time
 
    if (minorloop1(sinp,pinp))
        [fn, ln, sn, pn] = splitloop(sinp,pinp);    
        if(length(fn)~=0)
            minorloop{length(minorloop)+1} = fn;
            minorloop{length(minorloop)+1} = sn{1};
            minortime{length(minortime)+1} = ln;
            minortime{length(minortime)+1} = pn{1};
        end
    else
        if(length(sinp)~=0)
            minorloop{length(minorloop)+1} = sinp;
            minortime{length(minortime)+1} = pinp;
        end
    end
    
end

%***********************************************************************************************
%***********************************************************************************************


%=======================================================
%calculate loss for loop segment using improved equation
%=======================================================
function pseg = calcseg(t,B,alpha,beta,k1,a)

                            %a is fraction of Bpp used, 1-a is fraction of original
bma1a=(beta-alpha)*(1-a);   %exponent of B(t)
bmaa=(beta-alpha)*a;        %exponent of Bpp

Bpp=max(B)-min(B);

if sum(B<0)>0
	[t,B]=makepositive(t,B);
end

len=length(t);
T=t(len)-t(1);
deltaB=abs(B(2:len)-B(1:len-1));
deltat=(t(2:len)-t(1:len-1));
dBdt=deltaB./deltat;
m1=dBdt.^(alpha-1);
m2=abs((B(2:len)).^(bma1a+1)-(B(1:len-1)).^(bma1a+1));
pseg=k1/T/(bma1a+1)*sum(m1.*m2)*Bpp^bmaa;


%==========================================================================
%Convert a piecewise-linear waveform w (which must be a vector)
%represented by the points w at times t, to a piecewise-linear waveform with 
%points at any zero crossings.  Thus, any segment of the returned waveform is all 
%positive or all negative.  If w is a matrix, use makepositiveM.
%==========================================================================
function [t,w] = makepositive(t,w)
len = length(t);
cross = ( ( w(1:len-1).*w(2:len) ) < 0 );
loopnumber=len-1+sum(cross);
for i = 1:loopnumber
   if cross(i)
      len=length(w);
      tcross = w(i)*(t(i+1) - t(i))/(w(i)-w(i+1)) + t(i);
      t = [ t(1:i) tcross t(i+1:len)];
      w = [  w(1:i) 0 w(i+1:len)];
      cross = [ cross(1:i) 0 cross(i+1:len-1)];
   end
 end  
w=abs(w);



%==================
%minorloop function
%==================
function [ml] = minorloop1(s,p)
qr = length(s);
peak = -1;
prevslope = 0;
for i = 2 : qr
    slope = (s(1,i-1)-s(1,i))/(p(1,i-1)-p(1,i));
    if (slope <=0 & prevslope >=0)
        peak = peak + 1;
    end 
    if (prevslope <=0 & slope >=0)
        peak = peak + 1;
    end 
    prevslope = slope;
end 
if (peak > 2)
    ml = 1;
else
    ml = 0;
end