view release on metacpan or  search on metacpan

lib/ICC/Shared.pm  view on Meta::CPAN

			} else {
				
				# use lower source value
				$vector_out->[$i] = $vector_in->[$low];
				
			}
			
		}
		
	}

	# return
	return($vector_out);

}

# compute cubic spline interpolation matrix
# input/output range structure: [start_nm, end_nm, increment]
# optional extrapolation method is 'copy' or 'linear', none returns zeros
# parameters: (input_range, output_range, [extrapolation_method])
# returns: (interpolation_matrix)
sub cspline_matrix {

	# get parameters
	my ($range_in, $range_out, $ext) = @_;

	# local variables
	my ($ix, $ox, $rhs, $info, $derv, $mat);
	my ($w, $low, $t, $tc, $h00, $h01, $h10, $h11);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# compute input vector size from range
	$ix = round(($range_in->[1] - $range_in->[0])/$range_in->[2]);

	# verify input range
	($ix > 0 && abs($ix * $range_in->[2] - $range_in->[1] + $range_in->[0]) < 1E-12 && $range_in->[2] > 0) || croak('invalid input range');

	# compute output vector size from range
	$ox = round(($range_out->[1] - $range_out->[0])/$range_out->[2]);

	# verify output range
	($ox > 0 && abs($ox * $range_out->[2] - $range_out->[1] + $range_out->[0]) < 1E-12 && $range_out->[2] > 0) || croak('invalid output range');

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {
		
		# make rhs matrix (filled with zeros)
		$rhs = ICC::Support::Lapack::zeros($ix + 1);
		
		# for each row
		for my $i (1 .. $ix - 1) {
			
			# set diagonal values
			$rhs->[$i - 1][$i] = 3;
			$rhs->[$i + 1][$i] = -3;
			
		}
		
		# set endpoint values
		$rhs->[0][0] = -3;
		$rhs->[1][0] = -3;
		$rhs->[$ix][$ix] = 3;
		$rhs->[$ix - 1][$ix] = 3;
		
		# solve for derivative matrix
		($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
		
		# make matrix of zeros
		$mat = ICC::Support::Lapack::zeros($ox + 1, $ix + 1);
		
	# otherwise, use Math::Matrix package (slow)
	} else {
		
		# make rhs matrix (fill with zeros)
		$rhs = bless([map {[(0) x ($ix + 1)]} (0 .. $ix)], 'Math::Matrix');
		
		# for each row
		for my $i (1 .. $ix - 1) {
			
			# set diagonal values
			$rhs->[$i - 1][$i] = 3;
			$rhs->[$i + 1][$i] = -3;
			
		}
		
		# set endpoint values
		$rhs->[0][0] = -3;
		$rhs->[1][0] = -3;
		$rhs->[$ix][$ix] = 3;
		$rhs->[$ix - 1][$ix] = 3;
		
		# solve for derivative matrix
		$derv = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
		
		# make matrix of zeros
		$mat = [map {[(0) x ($ix + 1)]} (0 .. $ox)];
		
	}

	# for each output
	for my $i (0 .. $ox) {
		
		# compute wavelength
		$w = $range_out->[0] + $i * $range_out->[2];
		
		# if wavelength < start of source
		if ($w < $range_in->[0]) {
			
			# if extrapolation defined
			if (defined($ext)) {
				
				# if linear extrapolation
				if ($ext eq 'linear') {
					
					# compute ratio
					$t = ($w - $range_in->[0])/$range_in->[2];
					
					# for each input
					for my $j (0 .. $ix) {
						
						# set element to linear value
						$mat->[$i][$j] = ($j == 0 ? 1 : 0) + $derv->[0][$j] * $t;
						
					}
					
				# if copy extrapolation
				} elsif ('copy') {
					
					# set element to first knot
					$mat->[$i][0] = 1;
					
				} else {
					
					# error
					croak('invalid extrapolation type');
					
				}
				
			}
			
		# if wavelength > end of source
		} elsif ($w > $range_in->[1]) {
			
			# if extrapolation defined
			if (defined($ext)) {