function [XIB0,LINKS,P1,P2,REST,STIFF,dq,NIBS,HINGES,B1,B2,B3,BENDSTIFF,C,YT,TARGETSTIFF] = ...
                 ibInit(h,Slink,Sbend,Starget,demonumber);
%
% 
% Sets up initial IB points and data structures needed
% to evaluate IB forces
%
% inputs:  grid spacing h and spring stiffness Slink 
%
% outputs:  XIB0        initial IB point locations 
%           LINKS       list of link indices
%           P1,P2       indices IB points joined by link l
%           REST        rest length of link l
%           STIFF       stiffness of link l
%           dq          parameter increment
%           demonumber  choice of demo to run
%
%
% make provision for storing number of IB points in each of up to 10 IB objects
    NIBS = zeros(10,1);
%
%      
% DEMO 1 SETUP ************************************************************
%
    if demonumber == 1,

       disp(' demo = ellipse')
%
% ellipse  [a*cos(2*pi*theta)+0.5 b*sin(2*pi*theta)+0.5]
%
% compute arc length of ellipse
%
    theta = 2*pi*[0:.001:0.25];
    ltheta = length(theta)
    index   = [2:1:ltheta];
    indexm1 = [1:1:ltheta-1];
    a=0.4;b=0.15;
    tmp1 = a*(cos(theta(index))-cos(theta(indexm1)));
    tmp2 = b*(sin(theta(index))-sin(theta(indexm1)));
    size(tmp1)
    size(tmp2)
    darc = sqrt(tmp1.*tmp1 + tmp2.*tmp2)
    arclength = 4*sum(darc)
%
% divide arclength by (h/2) to determine number of IB points on 
% ellipse
    NIB1 = round(arclength/(h/2))
%
% define radius of circle of area same as ellipse and define
% spring rest lengths as circumference of circle divided by NIB1
%
    radius = sqrt(a*b)
    dtheta = 2*pi/NIB1;
    dq     = dtheta*radius;

% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
    REST   = (dq/4)*ones(NIB1,1);
%%%    REST   = dq*zeros(NIB1,1);
    Slink/(dq*dq)
    STIFF  = (Slink/(dq*dq))*ones(NIB1,1);
    clear index;
    index = [0:1:NIB1-1]';
    XIB0 = [a*cos(dtheta*index)+0.5 b*sin(dtheta*index)+0.5];
%%    XIB0 = [a*cos(dtheta*index)+0.2 b*sin(dtheta*index)+0.8];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
    LINKS = [1:NIB1]';
    P1 = [1:NIB1]';
    P2 = [2:NIB1 1]';
%
% set up lists of 'hinges' (3-pt forces)
%
    HINGES = [1:NIB1]';
    B1 = [NIB1 1:NIB1-1]';
    B2 = [1:NIB1]';
    B3 = [2:NIB1 1]';
    BENDSTIFF = Sbend*ones(NIB1,1);
%%%    C = -(dq*dq)*sin(2*pi/NIB1)*ones(NIB1,1);
    C = (dq*dq)*sin(2*pi/NIB1)*ones(NIB1,1);
    disp('XIB0 '),size(XIB0)
    disp('C    '),size(C)
    pause
    NIBS(1) = NIB1;
% set target point information (here there are no target points desired)
    YT = zeros(NIB1,2);
    TARGETSTIFF = zeros(NIB1,1);
%
% end of demo 1 setup
%
    elseif demonumber == 2,
%      
% DEMO 2 SETUP ************************************************************
%      
     disp(' demo = twocircles')
%
% radii and centers of two circles
%
     radius1 = 0.2;
     radius2 = 0.2;
     xcenter1 = 0.3;
     ycenter1 = 0.3;
     xcenter2 = 0.7;
     ycenter2 = 0.7;
%
% circle1  [radius1*cos(2*pi*theta)+xcenter1 radius1*sin(2*pi*theta)+ycenter1]
% circle2  [radius2*cos(2*pi*theta)+xcenter2 radius2*sin(2*pi*theta)+ycenter2]
%
% set up circle 1
%
    arclength = 2*pi*radius1
    NIB1 = round(arclength/(h/2))
    dtheta = 2*pi/NIB1;
    dq     = dtheta*radius1;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
    REST   = (dq/4)*ones(NIB1,1);
    Slink/(dq*dq)
    STIFF  = (Slink/(dq*dq))*ones(NIB1,1);
    clear index;
    index = [0:1:NIB1-1]';
    XIB0 = [radius1*cos(dtheta*index)+xcenter1 ...
            radius1*sin(dtheta*index)+ycenter1];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
    LINKS = [1:NIB1]';
    P1 = [1:NIB1]';
    P2 = [2:NIB1 1]';
%
% set up lists of 'hinges' (3-pt forces)
%
    HINGES = [1:NIB1]';
    B1 = [NIB1 1:NIB1-1]';
    B2 = [1:NIB1]';
    B3 = [2:NIB1 1]';
    BENDSTIFF = Sbend*ones(NIB1,1);
%%%%    C = -(dq*dq)*sin(2*pi/NIB1)*ones(NIB1,1);
    C = (dq*dq)*sin(2*pi/NIB1)*ones(NIB1,1);
%
% set up circle 2
%
    arclength = 2*pi*radius2
    NIB2 = round(arclength/(h/2))
    dtheta = 2*pi/NIB2;
    dq     = dtheta*radius1;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
% resting lengths for circle 2 links
    RESTN   = (dq/4)*ones(NIB2,1);
%%%    RESTN   = dq*zeros(NIB2,1);
% append resting lengths for circle 2 links to those for circle 1
    REST = [REST;RESTN];
    Slink/(dq*dq)
% stiffnesses for circle 2 links
    STIFFN  = (Slink/(dq*dq))*ones(NIB2,1);
% append stiffnesses for circle 2 links to those for circle 1
    STIFF = [STIFF;STIFFN];
    clear index;
    index = [0:1:NIB2-1]';
% define circle 2 IB point locations
    XIB0N = [radius2*cos(dtheta*index)+xcenter2 ...
             radius2*sin(dtheta*index)+ycenter2];
% append IB locations for circle 2 links to those for circle 1
    XIB0 = [XIB0;XIB0N];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
% list of links for circle 2 links
    LINKSN = [(NIB1+1):(NIB1+NIB2)]';
% append circle 2 list of links to those for circle 1
    LINKS  = [LINKS;LINKSN];
% indices of IB points joined by links of circle 2
    P1N = [(NIB1+1):(NIB1+NIB2)]';
    P2N = [(NIB1+2):(NIB1+NIB2) NIB1+1]';
% append circle 2 linked IB point indices to those for circle 1
    P1 = [P1;P1N];
    P2 = [P2;P2N];
%
% set up lists of 'hinges' (3-pt forces)
%
% list of hinges for circle 2
    HINGESN = [(NIB1+1):(NIB1+NIB2)]';
% append circle 2 list of hinges to those for circle 1
    HINGES = [HINGES;HINGESN];
% indices of circle 2 IB points in hinges
    B1N = [(NIB1+NIB2) (NIB1+1):(NIB1+NIB2-1)]';
    B2N = [(NIB1+1):(NIB1+NIB2)]';
    B3N = [(NIB1+2):(NIB1+NIB2) (NIB1+1)]';
% append circle 2 hinged IB point indices to those for circle 1
    B1 = [B1; B1N];
    B2 = [B2; B2N];
    B3 = [B3; B3N];
% define circle 2 bending stiffness
    BENDSTIFFN = Sbend*ones(NIB2,1);
% append circle 2 bending stiffnesses to those for circle 1
    BENDSTIFF = [BENDSTIFF;BENDSTIFFN];
%%%%    CN = -(dq*dq)*sin(2*pi/NIB2)*ones(NIB2,1);
% define circle 2 hinge C values
    CN = (dq*dq)*sin(2*pi/NIB2)*ones(NIB2,1);
% append circle 2 hinge C values to those for circle 1
    C = [C;CN];
% store number of IB points in each circle
    NIBS(1) = NIB1;
    NIBS(2) = NIB2;
% set up target point data (here no target points desired)
    YT = zeros(NIB1+NIB2,2);
    TARGETSTIFF = zeros(NIB1+NIB2,1);
%
% end of demo 2 setup
    elseif demonumber == 3,
%      
% DEMO 3 SETUP ************************************************************
%            
% perastalsis demo
%
     disp(' demo = perastalsis')
%
% locations of top and bottom IB walls
%
     ytopwall    = 0.75;
     ybottomwall = 0.25;
%
% set up bottom wall
%
    NIB1 = 1/(h/2);
    dq     = h/2;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
    REST   = (dq/2)*ones(NIB1,1);
    Slink/(dq*dq)
    STIFF  = (Slink/(dq*dq))*ones(NIB1,1);
    STIFF(NIB1)=0.0 % temporary fix
    clear index;
    index = [1:1:NIB1]';
    XIB0 = [h/2*index ybottomwall*ones(NIB1,1)];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
    LINKS = [1:NIB1]';
    P1 = [1:NIB1]';
    P2 = [2:NIB1 1]';
%
% set up lists of 'hinges' (3-pt forces)
%
    HINGES = [1:NIB1]';
    B1 = [NIB1 1:NIB1-1]';
    B2 = [1:NIB1]';
    B3 = [2:NIB1 1]';
    BENDSTIFF = zeros(NIB1,1);
    C = zeros(NIB1,1);
%
% set up top wall
%
    NIB2 = NIB1;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
    RESTN   = (dq/2)*ones(NIB1,1);
    REST    = [REST;RESTN];
    STIFFN  = (Slink/(dq*dq))*ones(NIB1,1);
    STIFFN(NIB1)=0.0 % temporary fix
    STIFF   = [STIFF;STIFFN]
    clear index;
    index = [1:1:NIB1]';
    XIB0N  = [h/2*index ytopwall*ones(NIB1,1)];
    XIB0  = [XIB0;XIB0N];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
    LINKSN = [NIB1+1:2*NIB1]';
    LINKS  = [LINKS;LINKSN];
    P1N = [NIB1+1:2*NIB1]';
    P2N = [NIB1+2:2*NIB1 NIB1+1]';
    P1 = [P1;P1N];
    P2 = [P2;P2N];
%
% set up lists of 'hinges' (3-pt forces)
%
    HINGESN = [NIB1+1:2*NIB1]';
    HINGES  = [HINGES;HINGESN];
    B1N = [2*NIB1 NIB1+1:2*NIB1-1]';
    B2N = [       NIB1+1:2*NIB1  ]';
    B3N = [       NIB1+2:2*NIB1 1]';
    B1  = [B1;B1N];
    B2  = [B2;B2N];
    B3  = [B3;B3N];
    BENDSTIFFN = zeros(NIB1,1);
    BENDSTIFF = [BENDSTIFF;BENDSTIFFN];
    CN = zeros(NIB1,1);
    C = [C;CN];
%
% set up initial target point locations to coincide
% with initial wall IB points
% 
    YT(:,1) = XIB0(:,1);
    YT(:,2) = XIB0(:,2);
    TARGETSTIFF = Starget*ones(2*NIB1,1);
%
% set up circle
%
%%%%%%%
     radius = 0.2;
     xcenter = 0.3;
     ycenter = 0.5;
%
% circle  [radius*cos(2*pi*theta)+xcenter radius*sin(2*pi*theta)+ycenter]
%
    arclength = 2*pi*radius
    NIB3 = round(arclength/(h/2))
    dtheta = 2*pi/NIB3;
    dq     = dtheta*radius;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
% resting lengths for circle 2 links
    RESTN   = (dq/4)*ones(NIB3,1);
%%%    RESTN   = dq*zeros(NIB3,1);
% append resting lengths for circle links to those for walls
    REST = [REST;RESTN];
    Slink/(dq*dq)
% stiffnesses for circle links
    STIFFN  = (Slink/(dq*dq))*ones(NIB3,1);
% append stiffnesses for circle links to those for walls
    STIFF = [STIFF;STIFFN];
    clear index;
    index = [0:1:NIB3-1]';
% define circle 2 IB point locations
    XIB0N = [radius*cos(dtheta*index)+xcenter ...
             radius*sin(dtheta*index)+ycenter];
% append IB locations for circle links to those for walls
    XIB0 = [XIB0;XIB0N];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
% list of links for circle links
    LINKSN = [(2*NIB1+1):(2*NIB1+NIB3)]';
% append circle list of links to those for walls
    LINKS  = [LINKS;LINKSN];
% indices of IB points joined by links of circle
    P1N = [(2*NIB1+1):(2*NIB1+NIB3)]';
    P2N = [(2*NIB1+2):(2*NIB1+NIB3) 2*NIB1+1]';
% append circle 2 linked IB point indices to those for walls
    P1 = [P1;P1N];
    P2 = [P2;P2N];
%
% set up lists of 'hinges' (3-pt forces)
%
% list of hinges for circle
    HINGESN = [(NIB1+1):(NIB1+NIB3)]';
% append circle list of hinges to those for walls
    HINGES = [HINGES;HINGESN];
% indices of circle IB points in hinges
    B1N = [(NIB1+NIB3) (NIB1+1):(NIB1+NIB3-1)]';
    B2N = [            (NIB1+1):(NIB1+NIB3)]';
    B3N = [(NIB1+2):(NIB1+NIB3) (NIB1+1)]';
% append circle hinged IB point indices to those for walls
    size(B1)
    size(B1N)
    B1 = [B1; B1N];
    B2 = [B2; B2N];
    B3 = [B3; B3N];
% define circle 2 bending stiffness
    BENDSTIFFN = Sbend*ones(NIB3,1);
% append circle 2 bending stiffnesses to those for circle 1
    BENDSTIFF = [BENDSTIFF;BENDSTIFFN];
%%%%    CN = -(dq*dq)*sin(2*pi/NIB3)*ones(NIB3,1);
% define circle 2 hinge C values
    CN = (dq*dq)*sin(2*pi/NIB3)*ones(NIB3,1);
% append circle 2 hinge C values to those for circle 1
    size(C)
    size(CN)
    C = [C;CN];
% set up target point data (here no target points desired)
    range = (2*NIB1+1):(2*NIB1+NIB3);
    YT(range,1:2) = 0;
    TARGETSTIFF(range,1) = 0;

%%%%%%S

    NIBS(1) = NIB1;
    NIBS(2) = NIB1;
    NIBS(3) = NIB3;
%
% end of demo 3 setup
%
    elseif demonumber == 4,
%      
% DEMO 4 SETUP ************************************************************
%                  
%
% stationary ellipse and moving circle with background force
%
     disp(' demo = circle in ellipse')

%
% ellipse  [a*cos(2*pi*theta)+0.5 b*sin(2*pi*theta)+0.5]
%
% compute arc length of ellipse
%
    theta = 2*pi*[0:.001:0.25];
    ltheta = length(theta)
    index   = [2:1:ltheta];
    indexm1 = [1:1:ltheta-1];
    a=0.45;b=0.35;
    tmp1 = a*(cos(theta(index))-cos(theta(indexm1)));
    tmp2 = b*(sin(theta(index))-sin(theta(indexm1)));
    size(tmp1)
    size(tmp2)
    darc = sqrt(tmp1.*tmp1 + tmp2.*tmp2)
    arclength = 4*sum(darc)
%
% divide arclength by (h/2) to determine number of IB points on 
% ellipse
    NIB1 = round(arclength/(h/2))
%
% define radius of circle of area same as ellipse and define
% spring rest lengths as circumference of circle divided by NIB1
%
    radius = sqrt(a*b)
    dtheta = 2*pi/NIB1;
    dq     = dtheta*radius;

% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
    REST   = (dq/4)*ones(NIB1,1);
    Slink/(dq*dq);
    STIFF  = (Slink/(dq*dq))*ones(NIB1,1);
    clear index;
    index = [0:1:NIB1-1]';
    XIB0 = [a*cos(dtheta*index)+0.5 b*sin(dtheta*index)+0.5];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
    LINKS = [1:NIB1]';
    P1 = [1:NIB1]';
    P2 = [2:NIB1 1]';
%
% set up lists of 'hinges' (3-pt forces)
%
    HINGES = [1:NIB1]';
    B1 = [NIB1 1:NIB1-1]';
    B2 = [1:NIB1]';
    B3 = [2:NIB1 1]';
%%%%    BENDSTIFF = Sbend*ones(NIB1,1);
    BENDSTIFF = 0*ones(NIB1,1);
    C = (dq*dq)*sin(2*pi/NIB1)*ones(NIB1,1);
    disp('XIB0 '),size(XIB0)
    disp('C    '),size(C)
    pause
    NIBS(1) = NIB1;
% set target point information 
    YT(:,1) = XIB0(:,1);
    YT(:,2) = XIB0(:,2);
    TARGETSTIFF = Starget*ones(NIB1,1);
%
%
% set up circle
%
%%%%%%%
     radius  = 0.15;
     xcenter = 0.3;
     ycenter = 0.5;
%
% circle  [radius*cos(2*pi*theta)+xcenter radius*sin(2*pi*theta)+ycenter]
%
    arclength = 2*pi*radius
    NIB2 = round(arclength/(h/2))
    dtheta = 2*pi/NIB2;
    dq     = dtheta*radius;
% rest length (delta s from notes -- since dq is length, r=#, ds=dq)
% resting lengths for circle 2 links
    RESTN   = (dq/4)*ones(NIB2,1);
% append resting lengths for circle links to those for walls
    REST = [REST;RESTN];
    Slink/(dq*dq)
% stiffnesses for circle links
    STIFFN  = (Slink/(dq*dq))*ones(NIB2,1);
% append stiffnesses for circle links to those for walls
    STIFF = [STIFF;STIFFN];
    clear index;
    index = [0:1:NIB2-1]';
% define circle 2 IB point locations
    XIB0N = [radius*cos(dtheta*index)+xcenter ...
             radius*sin(dtheta*index)+ycenter];
% append IB locations for circle links to those for walls
    XIB0 = [XIB0;XIB0N];
%
% set up list of 'links' (springs) -- link l joins points with
% indices P1(l) and P2(l)
%
% list of links for circle links
    LINKSN = [(NIB1+1):(NIB1+NIB2)]';
% append circle list of links to those for walls
    LINKS  = [LINKS;LINKSN];
% indices of IB points joined by links of circle
    P1N = [(NIB1+1):(NIB1+NIB2)]';
    P2N = [(NIB1+2):(NIB1+NIB2) NIB1+1]';
% append circle 2 linked IB point indices to those for ellipse
    P1 = [P1;P1N];
    P2 = [P2;P2N];
%
% set up lists of 'hinges' (3-pt forces)
%
% list of hinges for circle
    HINGESN = [(NIB1+1):(NIB1+NIB2)]';
% append circle list of hinges to those for walls
    HINGES = [HINGES;HINGESN];
% indices of circle IB points in hinges
    B1N = [(NIB1+NIB2) (NIB1+1):(NIB1+NIB2-1)]';
    B2N = [            (NIB1+1):(NIB1+NIB2)]';
    B3N = [(NIB1+2):(NIB1+NIB2) (NIB1+1)]';
% append circle hinged IB point indices to those for walls
    size(B1)
    size(B1N)
    B1 = [B1; B1N];
    B2 = [B2; B2N];
    B3 = [B3; B3N];
% define circle 2 bending stiffness
    BENDSTIFFN = Sbend*ones(NIB2,1);
% append circle 2 bending stiffnesses to those for circle 1
    BENDSTIFF = [BENDSTIFF;BENDSTIFFN];
% define circle 2 hinge C values
    CN = (dq*dq)*sin(2*pi/NIB2)*ones(NIB2,1);
% append circle 2 hinge C values to those for circle 1
    size(C)
    size(CN)
    C = [C;CN];
% set up target point data (here no target points desired)
    range = (NIB1+1):(NIB1+NIB2);
    YT(range,1:2) = 0;
    TARGETSTIFF(range,1) = 0;
%%%%%%S
    NIBS(1) = NIB1;
    NIBS(2) = NIB2;
%
% end of demo 4 setup
%
    else
       disp(' NO SUCH DEMO NUMBER')
    end  

    return;