# $Source: /usr/pkg/lib/perl5/site_perl/5.8.0/Text/CSV/Munge/Strain.pm $
# $Date: 2006-08-22 $

package Text::CSV::Munge::Strain;

use strict;
use warnings;
use Carp qw(carp croak);
use Math::Trig;
use Text::CSV::Munge qw( quote_neatly );

our ($VERSION) = '$Revision: 0.01 $' =~ m{ \$Revision: \s+ (\S+) }xm;

# Package vars defined at BEGIN for building RE's needed later down.
use vars qw ( @qualities );

BEGIN {

	# Standardize hash key naming
	sub enforce_key_rules {
		my $key = shift;
		$key =~ s{^\047*}{}m; # Strip left quote.
		$key =~ s{\047*$}{}m; # Strip right quote.
		$key =~ s{_}{ }gm;
		my @words = split q{ }, $key;
		for (@words) {
			$_ = ucfirst $_;
		}
		$key = join q{ }, @words;
		return $key;
	}

	@qualities = (
		'type', 'gage factor', 'transverse sensitivity', 'ohms',
		'orientation', 'location', 'array'
	);

    # Standaradize the case of allowed strain gage hash keys.
	for (@qualities) { $_ = enforce_key_rules($_) }
} # BEGIN


# Confuses editor syntax highlighting.
my $ASCII_POUND = "\043";
my $ASCII_SQUOT = "\047";

# Cheat to avoid illegal division by zero.
my $LSB = 0x1;


# Store the essential column keys of N rosettes.
sub set_col_key_sg {
	shift; # Lose the class.
	my $obj = shift; # Passed through from Text::CSV::Munge.pm
	my $col_key = shift;

	my $value = shift;
	my $reg_ex = join q{|}, @qualities;
	if ( $col_key =~ m/^($reg_ex)$/im) {
		$obj->set_col_key( $col_key, $value, @_ );
	} else {
		croak "Oops! column key '$col_key' does not match 'm/^($reg_ex)$/im'. \n";
	}
	return 'Pthhht! to Perl::Critic';
}

############
# ROSETTES #
############

# Append a subarray converted from 3 ue's to max/min ue & angle.
sub rosette_rect {
	shift; # Lose the class.
	my $obj = shift; # Passed through from Text::CSV::Munge.pm
	my ( $eu_1, $eu_2, $eu_3, $poisson ) = @_;
	my ( $max_p_eu, $min_p_eu, $theta );    # All arefs
	my $id = 'RR ' . join '-', ( $eu_1, $eu_2, $eu_3 );
	carp q{Oops! Sub rosette_rect expects 3 channel (column) IDs }
		. q{followed by and Poisson's ratio, not '}
		. join ( ', ', @_ ) . "'\n"
		unless scalar @_ == 4 && $_[3] < 1;
	for my $i ( 0 .. $#{ $obj->get_key('Column Data Arefs')->[0] } ) {
		( $max_p_eu->[$i], $min_p_eu->[$i], $theta->[$i] ) =
			rect_calc(
				$obj->get_key('Column Data Arefs')->[$eu_1][$i],
				$obj->get_key('Column Data Arefs')->[$eu_2][$i],
				$obj->get_key('Column Data Arefs')->[$eu_3][$i],
				$obj->get_key('Column Hrefs')->[$eu_1]->{'Transverse Sensitivity'},
				$obj->get_key('Column Hrefs')->[$eu_2]->{'Transverse Sensitivity'},
				$obj->get_key('Column Hrefs')->[$eu_3]->{'Transverse Sensitivity'},
				$poisson
			);
	}
	for ( 'Max', 'Min', 'Deg' ) {
	   push @{ $obj->get_key('Column Names') }, quote_neatly("$id $_");
	}
	push @{ $obj->get_key('Column Data Arefs') }, ( \@$max_p_eu, \@$min_p_eu, \@$theta );

	# Associate these rosette channels with their component gages.
	# and those gage channels with these rosette channels.
	my $ptr = $#{ $obj->get_key('Column Names') };
	$obj->sg_set_col_key('array',"$eu_1, $eu_2, $eu_3", $ptr-2, $ptr-1, $ptr);
	$obj->sg_set_col_key('array', (join q{, }, ($ptr-2,$ptr-1, $ptr)), $eu_1, $eu_2, $eu_3);

	return 'Pthhht! to Perl::Critic';
}

# Append a subarray converted from 3 ue's to max/min ue & angle.
sub rosette_delta {
	shift; # Lose the class.
	my $obj = shift; # Passed through from Text::CSV::Munge.pm
	my ( $eu_1, $eu_2, $eu_3, $poisson ) = @_;
	my ( $max_p_eu, $min_p_eu, $theta );    # All arefs
	my $id = 'DR ' . join '-', ( $eu_1, $eu_2, $eu_3 );
	die q{Oops! Sub resolve_rosette_rect expects 3 channel (column) IDs }
		. q{followed by and Poisson's ratio, not '}
		. join ( ', ', @_ ) . "'\n"
		unless scalar @_ == 4 && $_[3] < 1;
	for my $i ( 0 .. $#{ $obj->get_key('Column Data Arefs')->[0] } ) {
		( $max_p_eu->[$i], $min_p_eu->[$i], $theta->[$i] ) =
			delta_calc(
				$obj->get_key('Column Data Arefs')->[$eu_1]->[$i],
				$obj->get_key('Column Data Arefs')->[$eu_2]->[$i],
				$obj->get_key('Column Data Arefs')->[$eu_3]->[$i],
				$obj->get_key('Column Hrefs')->[$eu_1]->{'Transverse Sensitivity'},
				$obj->get_key('Column Hrefs')->[$eu_2]->{'Transverse Sensitivity'},
				$obj->get_key('Column Hrefs')->[$eu_3]->{'Transverse Sensitivity'},
				$poisson
			);
	}
	for ( 'Max', 'Min', 'Deg' ) {
		push @{ $obj->get_key('Column Names') }, quote_neatly("$id $_");
	}
	push @{ $obj->get_key('Column Data Arefs') }, ( \@$max_p_eu, \@$min_p_eu, \@$theta );

	# Associate these rosette channels with their component gages.
	# and those gage channels with these rosette channels.
	my $ptr = $#{ $obj->get_key('Column Names') };
	$obj->sg_set_col_key('array',"$eu_1, $eu_2, $eu_3", $ptr-2, $ptr-1, $ptr);
	$obj->sg_set_col_key('array', (join q{, }, ($ptr-2,$ptr-1, $ptr)), $eu_1, $eu_2, $eu_3);

	return 'Pthhht! to Perl::Critic';
}

# Calcuate any rosette's 3 ue values given the principal strains and angle.
# Use to generate simulated data for algorithm testing.
sub retro_rosette {
	shift; # Lose the class.
	my ( $max, $min, $theta, $g2_angle, $g3_angle, $poisson, $kt ) = @_;

	( $max, $min ) = retro_corr_k( $max, $min, $kt, $poisson );

	$theta    = deg2rad( -$theta );    # Undo the reversal
	$g2_angle = deg2rad($g2_angle);
	$g3_angle = deg2rad($g3_angle);

	# Formula per Measurments Group TN-515 Page 4 Formulae 2a, 2b, 2c.
	my $foo = ( $max + $min ) / 2;
	my $bar = ( $max - $min ) / 2;
	my @e123;
	push @e123, $foo + $bar * cos 2 * $theta;
	push @e123, $foo + $bar * cos 2 * ( $theta + $g2_angle );
	push @e123, $foo + $bar * cos 2 * ( $theta + $g3_angle );
	return @e123;    # Give back ue values for grids 1, 2 and 3.
}

# Return max and min principal strains plus angle of max in degrees.
sub rect_calc {
	my ( $e1, $e2, $e3 ) = rect_corr_k(@_);    # Correct transverse sens error.

	# Formula per Measurments Group TN-515 Page 4 Formula 3.
	my $foo = ( $e1 + $e3 ) / 2;
	my $bar = ( ( ( $e1 - $e2 )**2 + ( $e2 - $e3 )**2 ) )**0.5 / 2**0.5;

	# Formula per Measurments Group TN-515 Page 4 Formula 5 (not 4).
	# Angle expressed from Grid 1 to axis of principal strain!
	# Use of Least Significant Bit prevents zero-divide error.
	my $rad = atan( ( 2 * $e2 - $e1 - $e3 ) / ( ( $e1 - $e3 ) + $LSB ) ) / 2;
	return ( $foo + $bar, $foo - $bar, rad2deg($rad) );
}

# Return max/min principal strains and angle of max.
sub delta_calc {
	my ( $e1, $e2, $e3 ) = delta_corr_k(@_);    # Correct transverse sens error

	# Formula per Measurments Group TN-515 Page 5 Formula 6.
	my $foo = ( $e1 + $_[1] + $_[2] ) / 3;
	my $bar = 2**0.5/3 * (($e1 - $e2)**2 + ($e2 - $e3)**2 + ($e3 - $e1)**2)**0.5;

	# Formula per Measurments Group TN-515 Page 5 Formula 8 (not 7).
	# Angle expressed from Grid 1 to axis of principal strain!
	# Use of Least Significant Bit prevents zero-divide error.
	my $rad = atan(3**0.5 * ($e2 - $_[2]) / ((2 * $e1 - $e2 - $e3) + $LSB)) / 2;
	return ( $foo + $bar, $foo - $bar, rad2deg($rad) );
}

# Correct for transverse sensitivity ex-post-facto after resolving principal
# strains from uncorrected rosettes of any type.
# Formula per Vishay Micro-Measurements TN-509, page 6, formulae 16 & 17.
# Not used because more accurate with three K values versus one K value.
sub corr_k {
	my ( $ep, $eq, $k, $p ) = @_;
	my @cor;
	push @cor, ( 1 - $p * $k ) / ( 1 - $k**2 ) * ( $ep - $k * $eq );  # max prin
	push @cor, ( 1 - $p * $k ) / ( 1 - $k**2 ) * ( $eq - $k * $ep );  # min prin
	return @cor;
}

# Reverse correction for transverse sensitivity of any rosette.
# Use for generating correctable fake data.
sub retro_corr_k {
	my ( $ep, $eq, $k, $p ) = @_;
	my @uncor;
	my $x = ( 1 - $p * $k ) / ( 1 - $k**2 );

	# Reverse engineered by author using simultaneous linear equations.
	# Given P, Q, x and k, solve for p and q
	#   p - kq = P/x
	# -kp +  q = Q/x
	push @uncor, $ep / $x + $k * ( $eq + $k * $ep ) / $x / ( 1 - $k**2 );
	push @uncor, ( $eq + $k * $ep ) / $x / ( 1 - $k**2 );
	return @uncor;
}

# Correct for transverse sensitivity.
# Formula per Vishay Micro-Measurements TN-509, page 5, formulae 10-12.
sub rect_corr_k {
	my ( $e1, $e2, $e3, $k1, $k2, $k3, $p ) = @_;
	my @cor;
	push @cor, rect_corr_sub( $k1, $p ) * ( $e1 - $k1 * $e3 );
	push @cor, rect_corr_sub( $k2, $p ) * ( $e2 - $k2 * ( $e1 + $e3 - $e2 ) );
	push @cor, rect_corr_sub( $k3, $p ) * ( $e3 - $k3 * $e1 );
	return @cor;
}

# Correct for transverse sensitivity.
# Formula per Vishay Micro-Measurements TN-509, page 6, formulae 13-15.
sub delta_corr_k {
	my ( $e1, $e2, $e3, $k1, $k2, $k3, $p ) = @_;
	my @cor;
	push @cor, delta_corr_sub( $e1, $e2, $e3, $k1, $p );
	push @cor, delta_corr_sub( $e2, $e3, $e1, $k2, $p );
	push @cor, delta_corr_sub( $e3, $e1, $e2, $k3, $p );
	return @cor;
}

# Called by &rect_corr_k internally.
sub rect_corr_sub {
	my ( $k, $p ) = @_;
	return ( 1 - $p * $k ) / ( 1 - $k**2 );
}

# Called by &delta_corr_k internally.
sub delta_corr_sub {
	my ( $e1, $e2, $e3, $k, $p ) = @_;
	return ( 1 - $p * $k ) / ( 1 - $k**2 ) *
	  ( ( 1 + $k / 3 ) * $e1 - 2 / 3 * $k * ( $e2 + $e3 ) );
}

1;

__END__

=head1 NAME

Text::CSV::Munge::Strain.pm

=head1 VERSION

Version 0.01

=head1 SYNOPSIS

This module recieves its method calls as a pass-through from its parent module.
So your C<use> statment calls only that one, not this one. You also only refer
to that one, as below...

Establish new object.

C<use Text::CSV::Munge; >

C<my $strains = Text::CSV::Munge-E<gt>new();>

Let us assume we already have a CSV file containing seven columns and 500 lines.
Say that channel zero (1st column in the CSV) represents time as the X axis and
that channels 1 through 6 (columns 2 through 7 of the CSV) represent readings in
microstrain from two strain gage rosettes, one rectangular and one delta.

Assuming thus, we read in that file like so.

C<$strains-E<gt>merge_csvs( /some/path/sg_data.csv', []);>

We are assuming the input CSV is raw data. So let us add to that the usual
qualities expected to be known for each of our 6 gages. Since they are specific
to a given channel, and channels are represented in a CSV by single columns, we
assign these qualities by setting a key in a hash for that column, like so...

C<$strains-E<gt>set_col_key('ohms', 120, 1 .. 6);>

C<$strains-E<gt>set_col_key('gage factor', 2.01, 1 .. 6);>

C<$strains-E<gt>set_col_key('transverse sensitivity', 0.015, 1 .. 6);>

Return the established column keys of any channel (CSV column).

And then just to check we'll see how that took by asking the module to describe
those six columns by reading back a synopsis of all their column hash keys,
like so...

C<print $strains-E<gt>describe(1 .. 6);>

All six strain gages now have have the necessary column key values assigned
so that we may not proceed to calcuating their values into a combined, total
strain reading complete with angle in degrees relative to gage 1 of the trio.

Our first three strain gage channels comprise a rectangular rosette. So let
us now append three new channels (future CSV columns) describing the max and
min principal strains and the angle in degrees, like so...

C<$strains-E<gt>sg_rosette_rect( 1, 2, 3, $poisson );>

Thus are channels 7 through 8 (future CSV columns 8 through 11) now exist,
with channel names (future 1st-row column headings) appropriate to their data
type. Further, the three original strain channels have been associated to their
derivitive values by a new column hash key "array" whose value is a list of the
derived channel (0-based) numbers. The derived channels have a similar key
listing back to the source channel numbers.

Thus, with channels 7, 8 and 9 being derived channels, you may determine their
source channels like so...

C<my @source_chans = $strains-E<gt>get_col_key('array', 7);>

...or likewise going the other way, like so...

C<my @derived_chans = $strains-E<gt>get_col_key('array', 1);>

Let us now do likewise for our 2nd trio of strain gages for the delta rosette,
like so...

C<$strains-E<gt>sg_rosette_delta( 4, 5, 6, $poisson );>

...whereby we obtain similar results, except for its use of formula specific to
the differing geometry. So having called thus, three new channels are created
just as for the prior example.

We'll write out our calcuations saving only time and the two rosettes.

C<$strains->write_csv( /some/path/sg_and_rosette_data.csv', 0, 7 .. 12 );>

=head1 DESCRIPTION

Mostly for calculating the maximum and minimum principal strains from delta and
rectangular strain gage rosettes.

All formula were taken from online documentation by Vishay Measurements Group:

F<http://www.vishay.com/brands/measurements_group/guide/tn/tn515/515drnj.htm>

=over 4

=item Transverse Sensitivity

Calculations made for strain gage rosettes include correction for transverse
sensitivity using formulae garnered from Vishay Micro-Measurements Tech Note
TN-509 as detailed by page and formula number within related subroutines.

Rectangular and delta rosettes are treated separately, each with its own unique
correction formula particular to that type with separate GF and Kt for each
gage.

A generic retro-rosette subroutine generates uncorrected three-gage output
from corrected max/min principal strains. This method is of use mainly for
testing and debugging.

=item Wheatstone Bridge Nonlinearity

Not yet implementted. Will be per Vishay Micro-Measurments Tech Note
TN-507-1 when I get to it.

=back

=head1 BUGS AND LIMITATIONS

This program is free software. It carries absolutely no warranties or guarantees
of any kind (expressed, implied, or even vaguely hinted at).

=head1 AUTHOR

Gan Uesli Starling <F<gan@starling.us>>

=head1 LICENSE AND COPYRIGHT

Copyright (c) 2006 Gan Uesli Starling. All rights reserved.

This program is free software; you may redistribute and/or modify it under the
same terms as Perl itself.

=cut