use strict;
use lib '/people/sche155/Eye/create_sphere2/';
use generate_cmgui_mesh qw(generate_cmgui_mesh);
use Math::Trig qw(pi asin);


my $node_filename = $ARGV[0];
my $element_filename = $ARGV[1];

my $xi1_discretization = 8;
my $xi2_discretization = 8;
my $xi3_discretization = 8;

my $pi = pi;

my @centre = (307.6, 265.0, -119.2);

my @axis1 = (12.3, 0.0, 0.0);
my @axis2 = (0.0, 13.3, 0.0);
my @axis3 = (1.3, 0.0, 11.9);

#If $generate_cubics is non zero then the resulting mesh will be cubic Hermite, otherwise linear.
my $generate_cubics = 1;

my $group_name = "sclera";

sub location_function
{
  my ($xi1_step, $xi2_step, $xi3_step) = @_;

  my @coordinates;
  my @derivatives;

  my $xi3 = $xi3_step / $xi3_discretization;
  my $phi = asin(2.0 * $xi3 - 1.0);

  my $xi2 = $xi2_step / $xi2_discretization;
  my $psi = $pi * $xi2;

  my $xi1 = $xi1_step / $xi1_discretization;
  my $theta = $pi * $xi1;

  for (my $i = 0 ; $i < 3 ; $i++)
	 {
		$coordinates[$i] = $centre[$i] +
		  sin ($theta) * (2.0 * $xi2 - 1.0) * cos($phi) * $axis1[$i] +
			 (cos ($theta) * (1.0 - sin($psi)) + (1.0 - 2.0 * $xi1) * sin($psi)) * cos($phi) * $axis2[$i] +
				sin($phi) * $axis3[$i];

		# d(theta)/ds1
		my $dtheta = $pi / $xi1_discretization;

		# d(psi)/ds2
		my $dpsi = $pi / $xi2_discretization;

		# d(phi)/ds3
                my $dphi;
		if (($xi3 == 1.0) || ($xi3 == 0.0))	# at the top and bottom of sphere
		{
		   $dphi = 1.0;
		}
		else					# the rest of the sphere
		{
                   $dphi = 2.0 / sqrt (1.0 - (2.0 * $xi3 - 1.0) * (2.0 * $xi3 - 1.0)) / $xi3_discretization;
		}

		# d/ds1
		$derivatives[$i * 7] =
		  $dtheta * cos ($theta) * (2.0 * $xi2 - 1.0) * cos($phi) * $axis1[$i] +
			 (- $dtheta * sin ($theta) * (1.0 - sin($psi)) - 2.0 / $xi1_discretization * sin($psi)) * cos($phi) * $axis2[$i];

		# d/ds2
		$derivatives[$i * 7 + 1] =
		  2.0 / $xi2_discretization * sin ($theta) * cos($phi) * $axis1[$i] +
			 $dpsi * cos ($psi) * (- cos ($theta) + (1.0 - 2.0 * $xi1)) * cos($phi) * $axis2[$i];

		# d2/ds1ds2
		$derivatives[$i * 7 + 2] =
		  $dtheta * cos ($theta) * 2.0 / $xi2_discretization * cos($phi) * $axis1[$i] +
			- $dpsi * cos ($psi)  * (- $dtheta * sin ($theta)  - 2.0 / $xi1_discretization) * cos($phi) * $axis2[$i];

		# d/ds3
		$derivatives[$i * 7 + 3] =
		  - $dphi * sin ($theta) * (2.0 * $xi2 - 1.0) * sin($phi) * $axis1[$i] +
			 - $dphi * (cos ($theta) * (1.0 - sin($psi)) + (1.0 - 2.0 * $xi1) * sin($psi)) * sin($phi) * $axis2[$i] +
				$dphi * cos($phi) * $axis3[$i];

		# d2/ds1ds3
		$derivatives[$i * 7 + 4] =
		  - $dtheta * $dphi * cos ($theta) * (2.0 * $xi2 - 1.0) * sin($phi) * $axis1[$i] +
			 - $dphi * (- $dtheta * sin ($theta) * (1.0 - sin($psi)) - 2.0 / $xi1_discretization * sin($psi)) * sin($phi) * $axis2[$i];

		# d2/ds2ds3
		$derivatives[$i * 7 + 5] =
		  - 2.0 / $xi2_discretization * $dphi * sin ($theta) * sin($phi) * $axis1[$i] +
			 - $dpsi * $dphi * (- cos ($theta) * cos($psi) + (1.0 - 2.0 * $xi1) * cos($psi)) * sin($phi) * $axis2[$i];

		# d3/ds1ds2ds3
		$derivatives[$i * 7 + 6] =
		  - $dtheta * 2.0 / $xi2_discretization * $dphi * cos ($theta) * sin($phi) * $axis1[$i] +
			 - $dpsi * $dphi * ($dtheta * sin ($theta) * cos($psi) - 2.0 / $xi1_discretization * cos($psi)) * sin($phi) * $axis2[$i];

	 }

  return(\@coordinates, \@derivatives);
}


generate_cmgui_mesh( $node_filename, $element_filename,
     $xi1_discretization, $xi2_discretization, $xi3_discretization,
     \&location_function, $generate_cubics, $group_name );