clear;
theta_min = 0;
theta_max = pi;
theta_del = (theta_max-theta_min)/50;
theta = theta_min : theta_del : theta_max;

phi_min = 0;
phi_max = pi;
phi_del = (phi_max-phi_min)/50;
phi = phi_min : phi_del : phi_max;

gamma_min = 1000;

for th = 1 : size(theta, 2),
  for ph = 1 : size(phi, 2),
    w(3) = cos(theta(th));
    w(2) = sin(theta(th)) * sin(phi(ph));
    w(1) = sin(theta(th)) * cos(phi(ph));
    gamma = 128.661245 * eye(3);
    
    for j = 1:3,
      for l = 1:3,
	for i = 1:3,
	  for k = 1:3,
	    
	    if ( (i==j) & (k==l) )
	      if (i==k) 
		gamma (j, l) = gamma (j, l) - 17.031 * w(i) * w(k); 
	      else 
		gamma (j, l) = gamma (j, l) - 22.166 * w(i) * w(k); 
	      end;
	      break;
	    end;
	    
	    if ( ( (i==k) & (j==l) ) | ( (i==l) & (j==k) ) )
	      gamma (j, l) = gamma (j, l) + 182.536 * w(i) * w(k); 
	    end;
	  
	  end;
	end;
      end;
    end;

%     gamma
    a = min(eig(gamma));
    if (a < gamma_min) gamma_min = a; end;
    
  end;
end;

gamma_min
