function [f, df] = targetFunc(S, dos_treatment, simset, vsample)

    x = S(1); 
    y = S(2);

    % [dosV1, dosV2] = dosU(dos_treatment,vsample(1),vsample(2),vsample(3),vsample(4));
    % [dosU1, dosU2] = dosU(dos_treatment,optset.fu.u1(x,y),optset.fu.u2(x,y),optset.fu.u3(x,y),optset.fu.u4(x,y));

    if size(vsample,2) == 4         % when 4x set of voltages is known
        vsamples(:,:,1) = vsample(1); 
        vsamples(:,:,2) = vsample(2); 
        vsamples(:,:,3) = vsample(3); 
        vsamples(:,:,4) = vsample(4); 
        dosV = dosU(dos_treatment,vsamples);
    elseif size(vsample,2) == 2     % when 2x set of (dosx,dosy) is known
        dosV{1} = vsample(:,1); 
        dosV{2} = vsample(:,2); 
    end

    fu(:,:,1) = simset.fu.u1(x,y);
    fu(:,:,2) = simset.fu.u2(x,y);
    fu(:,:,3) = simset.fu.u3(x,y);
    fu(:,:,4) = simset.fu.u4(x,y);

    dos = dosU(dos_treatment,fu);

    u1 = simset.fu.u1; 
    u2 = simset.fu.u2; 
    u3 = simset.fu.u3; 
    u4 = simset.fu.u4; 

    du1dx = simset.dfu.du1dx; 
    du2dx = simset.dfu.du2dx; 
    du3dx = simset.dfu.du3dx; 
    du4dx = simset.dfu.du4dx; 

    du1dy = simset.dfu.du1dy; 
    du2dy = simset.dfu.du2dy; 
    du3dy = simset.dfu.du3dy; 
    du4dy = simset.dfu.du4dy; 

    f1 = dos{1} - dosV{1}; 
    f2 = dos{2} - dosV{2}; 
    f = f1^2 + f2^2 ; % target func

    switch dos_treatment

        case 'delta/sum (4 electrodes)' 

            sumphi = u1(x,y) + u4(x,y) + u2(x,y) + u3(x,y); 
            dif1   = dos{1} * sumphi; 
            dif2   = dos{2} * sumphi; 

            d_dif1_dx = du1dx(x,y) + du4dx(x,y) - du2dx(x,y) - du3dx(x,y);
            d_dif2_dx = du1dx(x,y) + du2dx(x,y) - du4dx(x,y) - du3dx(x,y);
            d_sum_dx  = du1dx(x,y) + du4dx(x,y) + du2dx(x,y) + du3dx(x,y);

            df1dx = ( d_dif1_dx * sumphi - dif1 * d_sum_dx ) / sumphi^2;
            df2dx = ( d_dif2_dx * sumphi - dif2 * d_sum_dx ) / sumphi^2;

            dfdx = 2 * (f1 * df1dx  +  f2 * df2dx);

            d_dif1_dy = du1dy(x,y) + du4dy(x,y) - du2dy(x,y) - du3dy(x,y);
            d_dif2_dy = du1dy(x,y) + du2dy(x,y) - du4dy(x,y) - du3dy(x,y);
            d_sum_dy  = du1dy(x,y) + du4dy(x,y) + du2dy(x,y) + du3dy(x,y);

            df1dy = ( d_dif1_dy * sumphi - dif1 * d_sum_dy ) / sumphi^2;
            df2dy = ( d_dif2_dy * sumphi - dif2 * d_sum_dy ) / sumphi^2;

            dfdy = 2 * (f1 * df1dy  +  f2 * df2dy);

            df = [dfdx dfdy]'; 

        case 'delta/sum (2 electrodes)' 

            df1dx = ( (du1dx(x,y) - du3dx(x,y)) * (u1(x,y) + u3(x,y)) - (u1(x,y) - u3(x,y)) * (du1dx(x,y) + du3dx(x,y)) ) / ((u1(x,y) + u3(x,y))^2); 
            df1dy = ( (du1dy(x,y) - du3dy(x,y)) * (u1(x,y) + u3(x,y)) - (u1(x,y) - u3(x,y)) * (du1dy(x,y) + du3dy(x,y)) ) / ((u1(x,y) + u3(x,y))^2); 

            df2dx = ( (du2dx(x,y) - du4dx(x,y)) * (u2(x,y) + u4(x,y)) - (u2(x,y) - u4(x,y)) * (du2dx(x,y) + du4dx(x,y)) ) / ((u2(x,y) + u4(x,y))^2); 
            df2dy = ( (du2dy(x,y) - du4dy(x,y)) * (u2(x,y) + u4(x,y)) - (u2(x,y) - u4(x,y)) * (du2dy(x,y) + du4dy(x,y)) ) / ((u2(x,y) + u4(x,y))^2); 

            dfdx = 2*f1*df1dx + 2*f2*df2dx;
            dfdy = 2*f1*df1dy + 2*f2*df2dy;

            df = [dfdx dfdy]';         

        case 'delta/sum (diagonal)'

            df1dx_13 = ( (u1(x,y)+u3(x,y))*(du1dx(x,y)-du3dx(x,y)) - (u1(x,y)-u3(x,y))*(du1dx(x,y)+du3dx(x,y)) ) / (u1(x,y)+u3(x,y))^2; 
            df1dx_24 = ( (u2(x,y)+u4(x,y))*(du2dx(x,y)-du4dx(x,y)) - (u2(x,y)-u4(x,y))*(du2dx(x,y)+du4dx(x,y)) ) / (u2(x,y)+u4(x,y))^2; 
            df1dy_13 = ( (u1(x,y)+u3(x,y))*(du1dy(x,y)-du3dy(x,y)) - (u1(x,y)-u3(x,y))*(du1dy(x,y)+du3dy(x,y)) ) / (u1(x,y)+u3(x,y))^2; 
            df1dy_24 = ( (u2(x,y)+u4(x,y))*(du2dy(x,y)-du4dy(x,y)) - (u2(x,y)-u4(x,y))*(du2dy(x,y)+du4dy(x,y)) ) / (u2(x,y)+u4(x,y))^2; 

            df2dx_13 = ( (u1(x,y)+u3(x,y))*(du1dx(x,y)-du3dx(x,y)) - (u1(x,y)-u3(x,y))*(du1dx(x,y)+du3dx(x,y)) ) / (u1(x,y)+u3(x,y))^2;
            df2dx_24 = ( (u2(x,y)+u4(x,y))*(du2dx(x,y)-du4dx(x,y)) - (u2(x,y)-u4(x,y))*(du2dx(x,y)+du4dx(x,y)) ) / (u2(x,y)+u4(x,y))^2; 
            df2dy_13 = ( (u1(x,y)+u3(x,y))*(du1dy(x,y)-du3dy(x,y)) - (u1(x,y)-u3(x,y))*(du1dy(x,y)+du3dy(x,y)) ) / (u1(x,y)+u3(x,y))^2; 
            df2dy_24 = ( (u2(x,y)+u4(x,y))*(du2dy(x,y)-du4dy(x,y)) - (u2(x,y)-u4(x,y))*(du2dy(x,y)+du4dy(x,y)) ) / (u2(x,y)+u4(x,y))^2; 

            dfdx = f1*(df1dx_13 - df1dx_24) + f2*(df2dx_13 + df2dx_24);
            dfdy = f1*(df1dy_13 - df1dy_24) + f2*(df2dy_13 + df2dy_24);

            df = [dfdx dfdy]'; 
            
        otherwise 
            disp(['Voltage inversion is not supported for ',dos_treatment,' DOS treatment.'])
            return
    end
end

