#######################################################################
#
# EllipticNomogramScaling.tcl
# Version: February 21, 2008
#
# Ron Doerfler
# doerfpub@myreckonings.com
# See http://www.myreckonings.com/wordpress/
#
# Calculates labels for the T,W,O scales on the upper and lower
# elliptic function curves y = +/- sqrt(4x^3-g2x-g3) for the equation
# relating the oxygen consumption of rainbow trout to the average
# weight and the water temperature:
# 
#                          O = k(T^n)(W^m)
#
#    where for T > 50 deg. F,
#        k = 3.05e-4
#        m = -0.138
#        n = 1.855
#
#    and for T < 50 deg. F,
#        k = 1.90e-6
#        m = -0.138
#        n = 3.130
#
# All input files are assumed to be in the same directory as this
# file. All output files are written to this same directory.
#
# The following files of generic labels of u,v and w for u+v+w=0 for
# g2=1200 and g3=-8500 are used as input:
#     positive_values_uppercurve_g2_1200_g3_-8500.txt (x,labelvalue)
#     negative_values_uppercurve_g2_1200_g3_-8500.txt (x,labelvalue)
#     positive_values_lowercurve_g2_1200_g3_-8500.txt (x,labelvalue)
#     negative_values_lowercurve_g2_1200_g3_-8500.txt (x,labelvalue)
#
# The following files are output for placing onto the elliptic
# function plot:
#     uppertickmarks.txt (xcoord,ycoord,angle) - angle uses scaled y
#     T_upper.txt (xcoord,ycoord,labelvalue,angle)
#     W_upper.txt (xcoord,ycoord,labelvalue,angle)
#     OX_upper_upper.txt (xcoord,ycoord,labelvalue,angle)
#     OX_upper_lower.txt (xcoord,ycoord,labelvalue,angle)
#     T_lower.txt (xcoord,ycoord,labelvalue,angle)
#     W_lower.txt (xcoord,ycoord,labelvalue,angle)
#     OX_lower_lower.txt (xcoord,ycoord,labelvalue,angle)
#     OX_lower_upper.txt (xcoord,ycoord,labelvalue,angle)
#
######################################################################

global env

# Define Weierstrass' Elliptic Function parameters.
set g2 1200
set g3 -8500

# Define maximum value of x to plot
set xmax 24

# Define scaling for x and y coordinates in typographic points (pt)
# NOTE that the TiKz plot command in LaTeX uses cm as units, so the
# x and y scaling for the plot command must be these values divided
# by the number of points in 1 cm = 28.45
set xscale 14
set yscale 0.70

# Define parameters for mapping T, O and W to u, v and w.
set upper_k 0.000305
set lower_k 0.00000190
set upper_n 1.855
set lower_n 3.130
set upper_m -0.138
set lower_m -0.138
set A 0.9725
set upper_u0 7.2478
set upper_w0 0.9725
set upper_v0 [expr -1.0 * ($upper_u0 + $upper_w0)]
set lower_u0 12.805
set lower_w0 0.335
set lower_v0 [expr -1.0 * ($lower_u0 + $lower_w0)]

# Define flags for indicating upper or lower curve
set upperflag 1
set lowerflag 0


proc calc_upper_y {xvalue} {
    # y = sqrt(4x^3-g2x-g3)
    global g2 g3
    set y [expr 4*($xvalue**3)-($g2*$xvalue)-$g3]
    if {$y < 0.0} {
        set y 0.0
    }
    return [string trimleft [format %.4f [expr sqrt($y)]] "0"]
}

proc calc_lower_y {xvalue} {
    # y = -sqrt(4x^3-g2x-g3)
    global g2 g3
    set y [expr 4*($xvalue**3)-($g2*$xvalue)-$g3]
    if {$y < 0.0} {
        set y 0.0
    }
    return [string trimleft [format %.4f [expr -1.0 * sqrt($y)]] "0"]
}

proc calc_upper_T {u} {
    # T = e^[A(u+u0)/n]
    global A upper_u0 upper_n
    return [string trimleft [format %#4.4g [expr exp($A*($u+$upper_u0)/$upper_n)]] "0"]
}

proc calc_lower_T {u} {
    # T = e^[A(u+u0)/n]
    global A lower_u0 lower_n
    return [string trimleft [format %#4.4g [expr exp($A*($u+$lower_u0)/$lower_n)]] "0"]
}

proc calc_upper_W {w} {
    # W = e^[A(w+w0)/m]
    global A upper_w0 upper_m
    return [string trimleft [format %#4.4g [expr exp($A*($w+$upper_w0)/$upper_m)]] "0"]
}

proc calc_lower_W {w} {
    # W = e^[A(w+w0)/m]
    global A lower_w0 lower_m
    return [string trimleft [format %#4.4g [expr exp($A*($w+$lower_w0)/$lower_m)]] "0"]
}

proc calc_upper_O {v} {
    # O = k*e^[-A(v+v0)]
    global A upper_k upper_v0
    return [string trimleft [format %#4.4g [expr $upper_k*exp(-1*$A*($v+$upper_v0))]] "0"]
}

proc calc_lower_O {v} {
    # O = k*e^[-A(v+v0)]
    global A lower_k lower_v0
    return [string trimleft [format %#4.4g [expr $lower_k*exp(-1*$A*($v+$lower_v0))]] "0"]
}

proc calc_upper_angle {xvalue} {
    # Tick/Label slope is perpendicular to y curve = -1/derivative = -1/(12x^2-g2)
    global xscale yscale g2
    if {$xvalue < -20} {
        return 180.0
    } elseif {$xvalue < -19.8} {
        return 170.0
    } elseif {$xvalue < -19.1} {
        return 160.0
    } else {
        set prev_x [expr $xvalue-1.0]
        set prev_y [calc_upper_y $prev_x]
        set prev_xscaled [expr $xscale * $prev_x]
        set prev_yscaled [expr $yscale * $prev_y]
        set next_x [expr $xvalue + 1.0]
        set next_y [calc_upper_y $next_x]
        set next_xscaled [expr $xscale * $next_x]
        set next_yscaled [expr $yscale * $next_y]
        set ave_slope [expr ($next_yscaled - $prev_yscaled)/($next_xscaled - $prev_xscaled)]
        return [format %4.4g  [expr 180 + (180.0/3.1416)*atan(-1.0/$ave_slope)]]
    }

}

proc get_angle {xvalue curve} {
    # Returns tick/label: curve = $upperflag or $lowerflag
    global upperflag lowerflag
    set infile2 [open [pwd]/uppertickmarks.txt r]
    while {[gets $infile2 angleline] >= 0} {
        set anglelist [split $angleline ","]
        set first [lindex $anglelist 0]
        if {$first == $xvalue} {
            set angle [lindex $anglelist 4]
            if {$curve == $lowerflag} {
                set angle [expr 180.0 - $angle]
            }
        }
    }
    close $infile2
    update
    update idletasks
    return $angle
}

proc inside_label_mods {curve} {
    # Special mods for inner labes in this SPECIFIC drawing to avoid
    # label collisions inside the curves when the separation narrows.
    # curve = $upperflag or $lowerflag
    global angle xcoord ycoord upperflag lowerflag
    if {($xcoord == 5) & ($curve == $upperflag)} {
        set angle [expr $angle - 10.0]
    } elseif {$xcoord == 6} {
        set angle [expr $angle - 10.0]
    } elseif {($xcoord == 7) & ($curve == $upperflag)} {
            set angle [expr $angle + 25.0]
    } elseif {($xcoord > 7.5) & ($xcoord < 12.5)} {
        set angle 90.0
        if {$curve == $upperflag} {
         set xcoord [expr $xcoord + 0.2]
        } else {
         set xcoord [expr $xcoord - 0.2]
        }
    } elseif {($xcoord == 13) & ($curve == $upperflag)} {
            set angle [expr $angle - 25.0]
    } elseif {$xcoord == 14} {
            set angle [expr $angle + 5.0]
    } elseif {$xcoord < -20} {
        if {$curve == $upperflag} {
         set ycoord [expr $ycoord + 4.0]
        } else {
         set ycoord [expr $ycoord - 4.0]
        }
    }
}


# TICKMARKS (x,y,angle)
# Populate file for upper tickmarks. The mark is perpendicular to the curve.
set infile [open [pwd]/positive_values_uppercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/uppertickmarks.txt}
set outfile [open [pwd]/uppertickmarks.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set tickangle [calc_upper_angle $xcoord] 
        if {$tickangle > 180.0} {
            set tickangle [expr $tickangle - 180.0]
        }
        if {$tickangle < 0.0} {
           set tickangle [expr $tickangle + 180.0]
        }
        set ycoord [calc_upper_y $xcoord]
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$tickangle"
    }
}
close $infile
close $outfile
update
update idletasks


# UPPER T-SCALE -- using positive upper curve values
set infile [open [pwd]/positive_values_uppercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/T_upper.txt}
set outfile [open [pwd]/T_upper.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,uvalue,Temp,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set uvalue [lindex $linelist 1]
        set ycoord [calc_upper_y $xcoord]
        set T [calc_upper_T $uvalue]
        set angle [get_angle $xcoord $upperflag]
        inside_label_mods $upperflag
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$uvalue,$T,$angle"
    }
}
close $infile
close $outfile
update
update idletasks

# UPPER W-SCALE -- using negative upper curve values
set infile [open [pwd]/negative_values_uppercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/W_upper.txt}
set outfile [open [pwd]/W_upper.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,wvalue,Weight,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set wvalue [lindex $linelist 1]
        set ycoord [calc_upper_y $xcoord]
        set W [calc_upper_W $wvalue]
        set angle [get_angle $xcoord $upperflag]
        inside_label_mods $upperflag
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]

        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$wvalue,$W,$angle"
    }
}
close $infile
close $outfile
update
update idletasks


# UPPER around to LOWER O-SCALE -- using positive upper curve
# values and positive lower curve values,
set infile [open [pwd]/positive_values_uppercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/OX_upper_upper.txt}
set outfile [open [pwd]/OX_upper_upper.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,vvalue,Oxygen,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set vvalue [lindex $linelist 1]
        set ycoord [calc_upper_y $xcoord]
        set OX [calc_upper_O $vvalue]
        set angle [get_angle $xcoord $upperflag]
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$vvalue,$OX,$angle"
    }
}
close $infile
close $outfile
update
update idletasks
set infile [open [pwd]/positive_values_lowercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/OX_upper_lower.txt}
set outfile [open [pwd]/OX_upper_lower.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,vvalue,Oxygen,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set vvalue [lindex $linelist 1]
        set ycoord [calc_lower_y $xcoord]
        set OX [calc_upper_O $vvalue]
        set angle [get_angle $xcoord $lowerflag]
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$vvalue,$OX,$angle"
    }
}
close $infile
close $outfile
update
update idletasks


# LOWER T-SCALE -- using negative lower curve values
set infile [open [pwd]/negative_values_lowercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/T_lower.txt}
set outfile [open [pwd]/T_lower.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,uvalue,Temp,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set uvalue [lindex $linelist 1]
        set ycoord [calc_lower_y $xcoord]
        set T [calc_lower_T $uvalue]
        # Adjust T=50 so it's not confused with upper T scale
        if {$T > 49.99} {
            set T 49.99
        }
        set angle [get_angle $xcoord $lowerflag]
        inside_label_mods $lowerflag
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$uvalue,$T,$angle"
    }
}
close $infile
close $outfile
update
update idletasks


# LOWER W-SCALE -- using negative lower curve values
set infile [open [pwd]/negative_values_lowercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/W_lower.txt}
set outfile [open [pwd]/W_lower.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,wvalue,Weight,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set wvalue [lindex $linelist 1]
        set ycoord [calc_lower_y $xcoord]
        set W [calc_lower_W $wvalue]
        set angle [get_angle $xcoord $lowerflag]
        inside_label_mods $lowerflag
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$wvalue,$W,$angle"
    }
}
close $infile
close $outfile
update
update idletasks


# LOWER around to UPPER OX-SCALE -- using positive lower curve
# values and positive upper curve values,
set infile [open [pwd]/positive_values_lowercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/OX_lower_lower.txt}
set outfile [open [pwd]/OX_lower_lower.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,vvalue,Oxygen,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set vvalue [lindex $linelist 1]
        set ycoord [calc_lower_y $xcoord]
        set OX [calc_lower_O $vvalue]
        set angle [get_angle $xcoord $lowerflag]
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$vvalue,$OX,$angle"
    }
}
close $infile
close $outfile
update
update idletasks
set infile [open [pwd]/positive_values_uppercurve_g2_1200_g3_-8500.txt r]
catch {file delete [pwd]/OX_lower_upper.txt}
set outfile [open [pwd]/OX_lower_upper.txt a]
puts $outfile "xcoord,ycoord,xscaled,yscaled,vvalue,Oxygen,angle"
while {[gets $infile pageline] >= 0} {
    set linelist [split $pageline ","]
    set xcoord [lindex $linelist 0]
    if {$xcoord <= $xmax} {
        set vvalue [lindex $linelist 1]
        set ycoord [calc_upper_y $xcoord]
        set OX [calc_lower_O $vvalue]
        set angle [get_angle $xcoord $upperflag]
        set xscaled [expr $xscale * $xcoord]
        set yscaled [expr $yscale * $ycoord]
        puts $outfile "$xcoord,$ycoord,$xscaled,$yscaled,$vvalue,$OX,$angle"
    }
}
close $infile
close $outfile
update
update idletasks

exit