function [f11_f, f21_f, f12_f, f22_f, f3_f, OD1_f, theta1_f, phi1_f, OD2_f, theta2_f, phi2_f] = nod_2O(xdata,Gamma,D_intra,D_iso,ydata)
%% Data Shaping
g=[xdata(:,1) xdata(:,2) xdata(:,3)];
G_mag= (xdata(:,4)./1000000); %(Tesla/micro meter)
Del=xdata(:,5);          %(in seconds)
del=xdata(:,6);          %(in seconds)
td=(Del-(del/3));        %(in seconds)
b=(td).*((Gamma.*del.*G_mag).^2); %(sec/mircometer^2)
[ydata1,sigma]=norm_meas_HCP(ydata,b);
 %% Estimation using GA
 %GA Module used from http://www.acse.dept.shef.ac.uk/cgi-bin/gatbx-download
NIND = 48;           
MAXGEN = 60;         % maximum Number of iterations
GGAP = 0.9;          % Generation gap, how many new individuals are created
NVAR = 8;            % Number of variables (OD theta phi v_ic) in each orientation
PRECI = 15;          % Precision of binary representation

% Build field descriptor
 FieldD =[15  15  15 15 15  15  15 15; 0.01  0.01  0.01 0.1 0.01  0.01  0.01 0.1;...
          0.98    pi    pi  1 0.98    pi    pi  1;1  1 1 1 1  1 1 1 ;0 0 0 0 0 0 0 0 ;1  1 1 1 1  1 1 1; 1   1 1 1 1   1 1 1];
% Initialise population
   Chrom = crtbp(NIND, NVAR*PRECI);

% Reset counters
   Best = NaN*ones(MAXGEN,1);	% best in current population
   gen = 0;                     % generational counter

% Evaluate initial population
par=bs2rv(Chrom,FieldD);
for i=1:NIND
Objv(i)=feval(@nod_eval_nlin_2O,par(i,:),D_intra,D_iso,G_mag,Del,del,g,b,sigma,ydata1);
end 
ObjV=Objv';
   Best(gen+1) = min(ObjV);
   while gen < MAXGEN,

    % Assign fitness-value to entire population
       FitnV = ranking(ObjV);

    % Select individuals for breeding
       SelCh = select('sus', Chrom, FitnV, GGAP);

    % Recombine selected individuals (crossover)
       SelCh = recombin('xovsp',SelCh,0.7);

    % Perform mutation on offspring
       SelCh = mut(SelCh);

    % Evaluate offspring, call objective function
    par1=bs2rv(SelCh,FieldD);
   
    for i=1:max(size(par1))
    ObjVsel(i)=feval(@nod_eval_nlin_2O,par1(i,:),D_intra,D_iso,G_mag,Del,del,g,b,sigma,ydata1); 
    end    
    ObjVSel = ObjVsel';  
    % Reinsert offspring into current population
       [Chrom ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);

    % Increment generational counter
       gen = gen+1;

    % Update display and record current best individual
       Best(gen+1) = min(ObjV);
   end 
[ObjVSel,ind]=sort(ObjVSel);
par1=par1(ind,:);
x1=par1(1,:);                % OD theta phi v_ic
%% Linear parameters fit using cvx
% toolbox available at http://cvxr.com/cvx/
[A11, A21, A12, A22, A3] = NOD_A_2O(x1,D_intra,D_iso,G_mag,Del,del,g,b);
A=[A11 A21 A12 A22 A3];
 cvx_begin quiet
     variable f11e 
     variable f21e 
     variable f12e 
     variable f22e 
     variable f3e 
     fe=[f11e f21e f12e f22e f3e]';
        minimize(norm((ydata1)- (A*fe)-sigma));   % for offset gaussian
       % minimize(norm((ydata1)- (A*fe)));   % for offset gaussian
     0.01 <= f11e <= 0.9
     0.01 <= f21e <= 0.9
     0.01 <= f12e <= 0.9
     0.01 <= f22e <= 0.9
     0.01 <= f3e <= 0.9
     f11e + f21e + f12e + f22e + f3e  == 1
  cvx_end
    f11ga=x1(4); f12ga=x1(8);
    x2=[((f11e+f11ga)/2) f21e ((f12e+f12ga)/2) f22e f3e x1(1) x1(2) x1(3) x1(5) x1(6) x1(7)];
%% Trust-region-reflective algorithm for non-linear data fitting
 opts = optimset('Display','off');
    lb=[0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01];
    ub=[0.9 0.9 0.9 0.9 0.9 0.99 pi pi 0.99 pi pi ];
    [xf,resnorm,~,exitflag,output] = lsqcurvefit(@noddi_eval_2O,x2,xdata,ydata1,lb,ub,opts);
f11_f=xf(1); f21_f= xf(2); f12_f=xf(3); f22_f= xf(4); f3_f=xf(5); OD1_f=xf(6); theta1_f=xf(7);  phi1_f=xf(8); OD2_f=xf(9); theta2_f=xf(10);  phi2_f=xf(11);

