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

# @ is used for arrays, $ is used for scalars
my $node_filename = "sclera_ori_layer.exnode";
my $element_filename = "sclera_ori_layer.exelem";

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

my $pi = pi;

# centre points of the eye ball
my @centre = (307.7, 265.0, -119.2);

# sclera thickness
my $thickness = 0.09;

# parameters for sclera layer
my @axis1 = (11.8, 0.0, 0.0);
my @axis2 = (0.0, 4.3, 10.9);
my @axis3 = (0.0, 11.9, -4.3);
#my @axis1 = (11.6, 0.0, 0.0);
#my @axis2 = (0.0, 12.6, 0.0);
#my @axis3 = (0.9, 0.0, 11.6);

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

my $group_name = "sclera";

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

# calculates the vector length of the sclera mesh
my $length1 = sqrt( $axis1[0] * $axis1[0] + $axis1[1] * $axis1[1] + $axis1[2] * $axis1[2] );
my $length2 = sqrt( $axis2[0] * $axis2[0] + $axis2[1] * $axis2[1] + $axis2[2] * $axis2[2] );
my $length3 = sqrt( $axis3[0] * $axis3[0] + $axis3[1] * $axis3[1] + $axis3[2] * $axis3[2] );

print "scleralength1 = $length1\n";
print "scleralength2 = $length2\n";
print "scleralength3 = $length3\n";


# new files for choroid layer
$node_filename = "choroid_ori_layer.exnode";
$element_filename = "choroid_ori_layer.exelem";

# calculate the new parameters for the choroid layer, the axis of the mesh
@axis1 = ( ($axis1[0] * (($length1 - $length1*$thickness) / $length1)) , ($axis1[1] * (($length1 - $length1*$thickness) / $length1)) , ($axis1[2] * (($length1 - $length1*$thickness) / $length1)) );

@axis2 = ( ($axis2[0] * (($length2 - $length2*$thickness) / $length2)) , ($axis2[1] * (($length2 - $length2*$thickness) / $length2)), ($axis2[2] * (($length2 - $length2*$thickness) / $length2)) );

@axis3 = ( ($axis3[0] * (($length3 - $length3*$thickness) / $length3)) , ($axis3[1] * (($length3 - $length3*$thickness) / $length3)), ($axis3[2] * (($length3 - $length3*$thickness) / $length3)) );

print "choroidaxis1 = @axis1\n";
print "choroidaxis2 = @axis2\n";
print "choroidaxis3 = @axis3\n";

# new thickness for choroid layer
$thickness = 0.0283;

$group_name = "choroid";

generate_cmgui_mesh( $node_filename, $element_filename,
     $xi1_discretization, $xi2_discretization, $xi3_discretization,
     \&location_function, $generate_cubics, $group_name );# @ is used for arrays, $ is used for scalars

# calculates the vector length of the choroid mesh
$length1 = sqrt( $axis1[0] * $axis1[0] + $axis1[1] * $axis1[1] + $axis1[2] * $axis1[2] );
$length2 = sqrt( $axis2[0] * $axis2[0] + $axis2[1] * $axis2[1] + $axis2[2] * $axis2[2] );
$length3 = sqrt( $axis3[0] * $axis3[0] + $axis3[1] * $axis3[1] + $axis3[2] * $axis3[2] );

print "clength1 = $length1\n";
print "clength2 = $length2\n";
print "clength3 = $length3\n";

# new files for the retina layer
$node_filename = "retina_ori_layer.exnode";
$element_filename = "retina_ori_layer.exelem";

# calculate the axis of the retina layer
@axis1 = ( ($axis1[0] * (($length1 - $length1*$thickness) / $length1)) , ($axis1[1] * (($length1 - $length1*$thickness) / $length1)) , ($axis1[2] * (($length1 - $length1*$thickness) / $length1)) );

@axis2 = ( ($axis2[0] * (($length2 - $length2*$thickness) / $length2)) , ($axis2[1] * (($length2 - $length2*$thickness) / $length2)), ($axis2[2] * (($length2 - $length2*$thickness) / $length2)) );

@axis3 = ( ($axis3[0] * (($length3 - $length3*$thickness) / $length3)) , ($axis3[1] * (($length3 - $length3*$thickness) / $length3)), ($axis3[2] * (($length3 - $length3*$thickness) / $length3)) );

print "retinaaxis1 = @axis1\n";
print "retinaaxis2 = @axis2\n";
print "retinaaxis3 = @axis3\n";

# new thickness for retina layer
$thickness = 0.029;

$group_name = "retina";

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


###############################################################################################################################
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 $xi2 = 1.0 - $thickness * $xi2;

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

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

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

		# d(psi)/ds2
		my $dxi2 =  - $thickness / $xi2_discretization;

		# d(phi)/ds3
                my $dphi;
		if (($xi3 == 1.0) || ($xi3 == 0.0))	# at the top and bottom of sphere
		{
		   $dphi = 0.75;
		}
		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) * $xi2 * cos($phi) * $axis1[$i]
			 - $dtheta * sin ($theta) * $xi2 * cos($phi) * $axis2[$i];

		# d/ds2
		$derivatives[$i * 7 + 1] =
		  $dxi2 * sin ($theta) * cos($phi) * $axis1[$i] +
			 $dxi2 * cos ($theta) * cos($phi) * $axis2[$i] +
				$dxi2 * sin($phi) * $axis3[$i];

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

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

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

		# d2/ds2ds3
		$derivatives[$i * 7 + 5] =
		  - $dxi2 * $dphi * sin ($theta) * sin($phi) * $axis1[$i]
		  	- $dxi2 * $dphi * cos ($theta) * sin($phi) * $axis2[$i] +
		  		$dxi2 * $dphi * cos($phi) * $axis3[$i];

		# d3/ds1ds2ds3
		$derivatives[$i * 7 + 6] =
		  - $dtheta * $dxi2 * $dphi * cos ($theta) * sin($phi) * $axis1[$i] +
			 $dtheta * $dxi2 * $dphi * sin ($theta) * sin($phi) * $axis2[$i];
	 }

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