
field_size = 180;

% do not have to be 180 units in Hz, 
% maybe some normalized units and then 
% scale it back to real f0 (Tilsen 2019's graph)?

nstep= 50;
begofactivation = 3;
endofactivation = 30;
timeforgest = 33;


t.targ = [150 80 80];

sim = Simulator();
sim.addElement(NeuralField('field f0', field_size, 10, -5, 4));
sim.addElement(LateralInteractions1D('f0 -> f0', field_size, 5, 10, 40,30, 0), ...
  'field f0', 'output', 'field f0', 'output');
% H+ L+ H-
sim.addElement(GaussStimulus1D('tone1',field_size, 5, 0, 10), ...
  [], [], 'field f0');
sim.addElement(GaussStimulus1D('tone2',field_size, 5, 0, 10), ...
  [], [], 'field f0');
sim.addElement(GaussStimulus1D('tone3',field_size, 50, 0, 10), ...
  [], [], 'field f0');


sim.tryInit();
sim.tryStep();

sim.init();

figure;

% only the planning course of the first H in HL or HH 
%is considered for now.

for istep = 1 : nstep
sim.step();
    act = sim.getComponent('field f0', 'activation');

    act_plus = act;
    act_plus (act_plus < 0) = 0;
    act = act - min(act);

    act_sum = sum(act);
    act_wgt = sum (act .* [1:field_size]);

    act_plus_sum = sum(act_plus);
    act_plus_wgt = sum (act_plus .* [1:field_size]);

    cen = act_wgt./act_sum;
    cen_plus = act_plus_wgt./act_plus_sum;
    centroid(istep)= cen_plus;
    

    if istep == begofactivation 
       sim.setElementParameters('tone1', 'amplitude', 10);
       sim.setElementParameters('tone1', 'position', t.targ(1));
       
       % H in HL
       sim.setElementParameters('tone2', 'amplitude',4);
       sim.setElementParameters('tone2', 'position', t.targ(2));
       sim.setElementParameters('tone3', 'amplitude', -8);
       sim.setElementParameters('tone3', 'position',90);
       
       %%% code for HH. no H-
       %sim.setElementParameters('tone2', 'amplitude',4);
       %sim.setElementParameters('tone2', 'position', t.targ(1));
       
       % The centroid shift is smaller than experimentally observed.
       % maybe the units of the field are not in Hz and add a scaling
       % factor to get real f0? Tilsen 2019.
    end

    
    if istep ==endofactivation 
        sim.setElementParameters('tone1', 'amplitude', 0);
        sim.setElementParameters('tone2', 'amplitude', 0);
        sim.setElementParameters('tone3', 'amplitude', 0);
    end 
    
        
    if istep == timeforgest
        plot(sim.getComponent('field f0', 'activation'), 'b');
        ylabel('act');
        if ~isnan(centroid(istep))
        xline(centroid(istep),'b'); %%centroid above 0
        end
        axis([0 180 -15 15])
        title({'Field activation', ['centroid:' num2str(centroid(istep))]})
        pause(0.2)
    end
        
end

%figure;
%plot (centroid, '-ob');
