function [x0 A b Aeq beq lb ub] = findQTlin1(s,freq)

% [x0 A b Aeq beq lb ub] = findQTlin(s,freq)
%
% Input Variables
%   s - 
%   freq -
%
% Output Variables
%   x0 -
%   A -
%   b -
%   Aeq -
%   beq - 
%   lb -
%   ub- 
%
% Description
%   Generates initial value and constraints for findQT.

% Created
%   4/22/2006, Richard J. Povinelli, Marquette University
%
% Modified
%   Mohamed A Mneimneh, Marquette University
% 
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).


s = s(:)'; %make sure s is a row vector

%right venticular endocardial excitation (RVen)
%right venticular epicardial excitation (RVep)
%left venticular endocardial excitation (LVen)
%left venticular epicardial excitation (LVep)
[regions FEATURES_PER_MODEL] = getRegions(); %heart regions modeled

LENGTH = length(regions) * FEATURES_PER_MODEL;

x0 = zeros(1, LENGTH);
lb = x0;
ub = x0;
[temp maxDeviationIndex] = max(abs(s));

ii2 = rand() * maxDeviationIndex; %initial random inflection
delta1 = rand() * freq * 0.1; %0.1s 
delta2 = rand() * freq * 0.5; %0.5s
ii1 = ii2 - delta2; %initial random inflection
x0(1:2) = [ii1 (ii1 + delta2)]; %inflection points for RVen
x0(3:4) = [(ii1 + delta1) (ii1 + delta2)]; %inflection points for RVep
x0(5:6) = [ii2 (ii2 + delta2)]; %inflection points for LVen
x0(7:8) = [(ii2 + delta1) (ii2 + delta2)]; %inflection points for LVep
x0(9:10) = [(ii2 + delta1) (ii2 + delta2)]; %inflection points for LVep

% x0(1:2) = [ii1 (ii1 + delta2)]; %inflection points for RVen
% x0(3:4) = [(ii1 + delta1) (ii1 + delta2)]; %inflection points for RVep
% x0(5:6) = [ii2 (ii2 + delta2)]; %inflection points for LVen
% x0(7:8) = [(ii2 + delta1) (ii2 + delta2)]; %inflection points for LVep

AV(1,:) = [1 -1 0 0 0 0 0 0 0 0]; bV(1) = 0;% RVen.c1 <= RVen.c2
AV(2,:) = [0 0 1 -1 0 0 0 0 0 0]; bV(2) = 0;% RVep.c1 <= RVep.c2
AV(3,:) = [0 0 0 0 1 -1 0 0 0 0]; bV(3) = 0;% RVep.c1 <= RVep.c2
AV(4,:) = [0 0 0 0 0 0 1 -1 0 0]; bV(4) = 0;% RVep.c1 <= RVep.c2
AV(5,:) = [0 0 0 0 0 0 0 0 1 -1]; bV(5) = 0;% RVep.c1 <= RVep.c2

AV(6,:) = [1 0 -1 0 0 0 0 0 0 0]; bV(6) = 0;% RVen.c1 <= RVep.c1
AV(7,:) = [0 1 -1 0 0 0 0 0 0 0]; bV(7) = 0;% RVep.c2 <= RVen.c2
AV(8,:) = [0 -1 1 0 0 0 0 0 0 0]; bV(8) = 0;% RVep.c2 <= RVen.c2

AV(9,:) = [0 0 1 0 -1 0 0 0 0 0]; bV(9) = 0;% RVen.c1 <= RVep.c1
AV(10,:) = [0 0 0 -1 0 1 0 0 0 0]; bV(10) = 0;% RVep.c2 <= RVen.c2

%limit time between inflections points to 50 - 300ms?

AV(11,:) = [0 0 0 0 0 0 1 0 -1 0]; bV(11) = 0;% LVen.c1 <= LVep.c1
AV(12,:) = [0 0 0 0 0 0 0 -1 0 1]; bV(12) = 0;% LVep.c2 <= LVen.c2

AV(13,:) = [0 0 0 0 1 0 -1 0 0 0]; bV(13) = 0;% RVep.c1 <= LVen.c1

AV(14,:) = [0 0 0 0 0 0 1 0 -1 0]; bV(14) = 0;% RVep.c1 <= LVen.c1
AV(15,:) = [0 0 0 0 0 0 0 -1 0 1]; bV(15) = 0;% LVep.c2 <= LVen.c2

lb(1:10) = zeros(10,1); %inflection point bounds
ub(1:2:10) = maxDeviationIndex * ones(5,1);
ub(2:2:10) = length(s) * ones(5,1);


x0(11:12) = [20 2]; %slope paramters for RVen
x0(13:14) = [20 6]; %slope paramters for RVep
x0(15:16) = [20 2]; %slope paramters for LVen
x0(17:18) = [20 6]; %slope paramters for LVep
x0(19:20) = [20 2]; %slope paramters for LVep
lb(11:20) = 0 * ones(10,1); ub(11:20) = 50 * ones(10,1); %slope bounds


x0(21) = 0.5* (max(s) - min(s)); %Magnitude RVen
x0(22) = 0.2* (max(s) - min(s)); %Magnitude RVep
x0(23) = max(s) - min(s); %Magnitude LVen
x0(24) = 0.2* (max(s) - min(s)); %Magnitude LVep
x0(25) = 0.2* (max(s) - min(s)); %Magnitude LVep
lb([21 22 24 25]) = zeros(4,1); ub([21 22 24 25]) = 2 * abs(max(s)) * ones(4,1); %mag bounds
lb(23) = 0.9* abs(max(s) - min(s)); ub(23) = 2 * abs(max(s) - min(s)); %LVen mag bounds



x0(24:25) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for RVen
x0(26:27) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for RVep
x0(28:29) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVen
x0(30:31) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVep
x0(32:33) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for RVen
x0(34:35) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for RVep
x0(36:37) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVen
x0(38:39) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVep
x0(40:41) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVen
x0(42:43) = [(rand() * freq * 0.2) (rand() * freq * 0.2)]; %Delays to positive and negative leads for LVep
lb(24:43) = zeros(20,1); ub(24:43) = freq * 0.25 * ones(20,1); %delay bounds




%form the big A matrix and b vector
A = blkdiag(AV, zeros(1,LENGTH-min(size(AV))));
b = [bV 0];
Aeq = [];
beq = [];