function [jav,divB] = curlB(Ba,Bb,Bc,Bd,ra,rb,rc,rd,y0)
%**************************************************************************
% CLUSTER project: CurlB                           Autor: Christian Maszl 
% Formulas: y0*Jav*(rab x rac) = dBab*rac-dBac*rab                        
%           div B*|rab*(rac x rad)| = |Sum(cyclic) dBab*(rac x rad)|    
%           a,b,c,d - 1,2,3,4                                            
% Unit System = [SI]                  Mail: e0025754@student.tuwien.ac.at 
%**************************************************************************
%
% Version: 1.0 - 20.08.2004
%
% This function computes the current density jav and divB from
% magnetometer data.

% Lock function in memory
mlock;
% Calculate average of B and convert units to SI
    dBab = (Ba-Bb)*1e-9; %[B]=T
    dBac = (Ba-Bc)*1e-9;
    dBad = (Ba-Bd)*1e-9;
    dBcb = (Bc-Bb)*1e-9;
    dBcd = (Bc-Bd)*1e-9;
% Calculate the differences and convert to SI
    rab  = (rb-ra)*1e+3; %[r]=m
    rac  = (rc-ra)*1e+3;
    rad  = (rd-ra)*1e+3;
    rcb  = (rb-rc)*1e+3;
    rcd  = (rd-rc)*1e+3;
    clear Ba Bb Bc Bd ra rb rc rd;
% Calculate the matrix components of f(B) 
    fB1 = dot(dBab,rac,2)-dot(dBac,rab,2); 
    fB2 = dot(dBac,rad,2)-dot(dBad,rac,2); 
    fB3 = dot(dBad,rab,2)-dot(dBab,rad,2);
    fB4 = dot(dBcb,rcd,2)-dot(dBcd,rcb,2);
% Calculate the matrix components of A (area)    
    A1  = -cross(rab,rac,2)*y0; 
    A2  = -cross(rac,rad,2)*y0; 
    A3  = -cross(rad,rab,2)*y0;
    A4  = -cross(rcb,rcd,2)*y0;
% Calculate the matrix components of V (volume)
    V   = abs(dot(rab,A2,2));
% Calculate the matrix components of fBV
    fBV=abs(dot(dBab,A2,2)+dot(dBad,A1,2)+dot(dBac,A3,2));
% Clear dispensable variables
    clear dBab dBac dBad dBcd dBcb;
    clear rab rac rad rcd rcb;
% Calculate <div(B)> and divide by y0 to get the same dimension
% as curlB
    divB=fBV./(V*y0);
% Preallocate the array jav
    jav =zeros(size(A1,1),3);
% Calculate the current density
for i=1:size(A1,1)
% Solve the overdetermined system of equations
jav(i,:) =([ A1(i,:);  A2(i,:);  A3(i,:); A4(i,:)]\ ...
           [fB1(i,:); fB2(i,:); fB3(i,:);fB4(i,:)])';
end
