#!/usr/bin/perl 
#
#
# Programm to calculate the magnetic field of conducting lines (in 2D)
#
# written by Michael Panhorst (Michael@Panhorst.com)

package EM;

our $VERSION = 0.25;

use strict;
use warnings;

use PDL ;
use PDL::Graphics::PGPLOT ;
use PDL::Graphics::LUT ;

use PDL::IO::Dumper ;



our @AllLines ;



##################################################################
# Check for Options and execute subroutines for the given switch

if (!@ARGV) { showhelp();}
elsif ($ARGV[0] eq "-h") { showhelp();}
elsif ($ARGV[0] eq "-f") { shift @ARGV; text_input($ARGV[0]) ;}
elsif ($ARGV[0] eq "-p") { shift @ARGV; image_input($ARGV[0]) ;}
else { print_error("unknown options: @ARGV"); }




##################################################################
# Some miscellaneous subroutines
################################

sub showhelp {
    print "\n" ;
    print "$0: Program to calculate the magnetic field of conducting lines (in 2D)\n" ;
    print "* * *   Version $VERSION, written by Michael Panhorst <Michael\@Panhorst.com>   * * *\n
" ;
    print "\n" ;
    print "    USAGE:\n" ;
    print "    ------\n" ;
    print "    $0 -f <Input text file>\n" ;
    print "    $0 -p <Input picture>\n\n" ;
    exit 1;
} ;


sub print_error {
    print "\nERROR: @_ \n" ;
    print "-----> showing help:\n" ;
    showhelp () ;
} ;


sub check_for_data {
    my $InputFile = shift ;

    print "\nChecking for actual data... \n" ;

    my $DataFile = $InputFile ;
    $DataFile =~ s/\..*$/.pdl/ ;


    if ( -r $DataFile ) {
	if ( -M $DataFile <= -M $InputFile ) {

	    print "-> FOUND! Displaying Data directly!\n\n" ;
	    my $Matrix = frestore($DataFile) ;

	    output($Matrix,$InputFile) ;

	}
    }

    print "-> No actual data found. Start calculating...\n\n" ;    
}





##################################################################
# additional calculating subroutines:


# Returns the vector length of (x,y) in the grid
sub grid_vec_length {
    my ( $x , $y ) = @_ ;
    return sqrt( $x * $x + $y * $y ) ; 
}


# Returns the real vector length of (x,y)
sub real_vec_length {
    my ( $x , $y ) = @_ ;
    return sqrt( $x * $x + $y * $y ) * $EM::GridWidth * 0.000000001 ; 
}

# Returns the vector direction of (x,y), divides by the length
sub vec_direction {
    my ( $x1 , $y1 , $x2 , $y2 , $width ) = @_ ;

    my $a = $x2 - $x1 ;
    my $b = $y2 - $y1 ;

    my $NIteration = grid_vec_length( $a , $b ) * $EM::Iteration;

    print_error("Vector with 0 length defined!\n") if $NIteration == 0 ;

    $a = $a / $NIteration ;
    $b = $b / $NIteration ;

    print_error("Width ($width) must be an integer!\n") if $width != int($width) ;

#    my $NoOfLines = $width * $EM::Iteration if $width ;
    my $NoOfLines = $width if $width ;

    return $a , $b , round( $NIteration ) , $NoOfLines ;
}


sub round {
    my($number) = shift;
    return int($number + .5 * ($number <=> 0));
}





sub calc_biot_savart {
    my $Distance = shift ;

    my $BSDistance = 2 * 0.0000001 * ( $EM::Current * 0.001 ) / ($Distance * 0.000000001) ;



#    our $BSThickness = 4 * 0.0000001 * ( $EM::Current * 0.001 ) / ($EM::Thickness * 0.000000001) ;
#    print "Biot-Savart for half thickness = $BSThickness [N/Am]\n" ;
    my $BSGrid = 2 * 0.0000001 * ( $EM::Current * 0.001 ) / ($EM::GridWidth * 0.000000001) ;
    print "Biot-Savart for Gridwidth = $BSGrid [N/Am]\n" ;



    return $BSDistance ;
}





##################################################################
# Main processing subroutines:
##############################


##################################################################
# TextFile-Input:


sub text_input {
    my $TextFile = shift ;

    print_error("Cannot read TextFile: $TextFile \n") unless -r $TextFile ;

    do $TextFile ;

    print "\nMatrix-size (X x Y): $EM::MatrixWidth x $EM::MatrixHeight \n";
    print "Current    = $EM::Current mA\n" ;
    print "mue0/4pi   = 10^-7 N/A^2\n" ;
    print "Grid width = $EM::GridWidth nm\n" ;
    print "$EM::Iteration Iterations per grid width\n" ;


    check_for_data($TextFile) ;



    # Fill the Matrix with zeroes:

    my $Matrix = zeroes ( $EM::MatrixHeight , $EM::MatrixWidth );


    foreach my $ConductingLines ( sort keys %$EM::Lines ) {

	my @StartPoint = ( $$EM::Lines{$ConductingLines}[0] , $$EM::Lines{$ConductingLines}[1] );
	my @dl = vec_direction (@{ $$EM::Lines{$ConductingLines} }) ;

	my @REALdl = @dl ;
	$_ = $_ * $EM::GridWidth * 0.000000001 foreach @REALdl ;


	print "Calculating: $ConductingLines => [ @{ $$EM::Lines{$ConductingLines} } ] with ", $dl[2], " * ", $dl[3]+1, "  lineparts\n" ;


	my $Konst = $EM::Current / 1000 * 0.0000001 ;


	$StartPoint[0] = $StartPoint[0] + $dl[1] * $EM::Iteration * $dl[3] / 2  if $dl[3] > 0 ;
	$StartPoint[1] = $StartPoint[1] - $dl[0] * $EM::Iteration * $dl[3] / 2  if $dl[3] > 0 ;

	my $PartCurrent = $EM::Current / ($dl[3]+1) if $dl[3] > 0 ;

	# Konstant = Current / 1000 (-> A) * 0.0000001 (mue0/4pi)
	$Konst = $PartCurrent / 1000 * 0.0000001 if $dl[3] > 0 ;

	for ( my $L = 0 ; $L <= $dl[3] ; $L++ ) {

	    my @BasePoint = @StartPoint ;		


	    print "StartPoint: @StartPoint \n";
#	    print "DL: @dl \n";
#	    print "Real-DL: @REALdl \n";


	    for ( my $n = 0 ; $n < $dl[2] ; $n++ ) {


#	    print "Part " , $n+1, " of ", $dl[2], " lineparts.\n" ;

		for ( my $Y = 0 ; $Y < $EM::MatrixHeight ; $Y++ ) {
		    
		    for ( my $X = 0 ; $X < $EM::MatrixWidth ; $X++ ) {

#		    print "X = $X \nY = $Y \n" ;

			my $r1 = $X - $BasePoint[0] ;
			my $r2 = $Y - $BasePoint[1] ;

			
			my $REALr1 = $r1 * $EM::GridWidth * 0.000000001 ;
			my $REALr2 = $r2 * $EM::GridWidth * 0.000000001 ;

			my $rlength = real_vec_length ( $r1, $r2 ) ;
			next unless $rlength ;

			my $Bz = $Konst * ( $REALdl[0] * $REALr2 - $REALdl[1] * $REALr1 ) / ($rlength * $rlength * $rlength) ;
			my $NewVal = at($Matrix, $Y, $X) + $Bz ;

			set ($Matrix, $Y, $X, $NewVal);

                # Check a special point (testing purpose only)
 		#    if ( $X == 2 && $Y == 2 ) {
	 	#	print "Basepoints: @BasePoint\n";
		#	print "dl: @dl \n" ;
		#	print "X = $X, Y = $Y \n" ;
		#	print "Bz = $Bz \n" ;
		#	print "rlength = $rlength \n\n" ;
		#    }


		    }
		
		}

		push @AllLines, round($BasePoint[0]) ;
		push @AllLines, round($BasePoint[1]) ;

#	    print "\n\n @BasePoint \n" ;
		$BasePoint[0] = $BasePoint[0] + $dl[0] ;
		$BasePoint[1] = $BasePoint[1] + $dl[1] ;

	    }
	   

 	    $StartPoint[0] = $StartPoint[0] - $dl[1] * $EM::Iteration ;
	    $StartPoint[1] = $StartPoint[1] + $dl[0] * $EM::Iteration ;
	}

    }

#    my $Check = at($Matrix, 2 , 4 ) ;
#    print "At(2,4) = $Check\n" ;


    my $DataFile = $TextFile ;
    $DataFile =~ s/\..*$/.pdl/ ;

    fdump($Matrix,$DataFile) ;


    output($Matrix,$TextFile) ;



}







##################################################################
# image Input (not yet implemented)



sub image_input {

    print "The -p flag is not yet implemented!\n";
    exit 1;

}









##################################################################
# What to do with the given Matrix (Output):
############################################


sub output {
    my $Matrix = shift ;
    my $InputFile = shift ;

    normalize_imag($Matrix,$InputFile);



#    print $Matrix ;
#    imag ($Matrix) ;
#    hold ;
#    cont($Matrix);



    exit 0 ;


}


##################################################################
# normalize your own way ...


sub normalize_imag{
    my $Matrix = shift ;
    my $InputFile = shift ;

    print $Matrix if $EM::MatrixHeight * $EM::MatrixWidth < 500 ;

    my $Min = 0;
    my $Max = 0 ;
    my $Mittelwert = 0 ;

    for ( my $Y = 0 ; $Y < $EM::MatrixHeight ; $Y++ ) {

	for ( my $X = 0 ; $X < $EM::MatrixWidth ; $X++ ) {
	    
#	    print "X = $X \nY = $Y \n" ;

	    $Max = at($Matrix, $Y, $X) if at($Matrix, $Y, $X) > $Max ;
	    $Min = at($Matrix, $Y, $X) if at($Matrix, $Y, $X) < $Min ;

	    $Mittelwert = $Mittelwert + abs(at($Matrix, $Y, $X)) ;



	}
    }


    ##################################################################
    # print values of one line
    my $OneX = 30 ; # middle
    my $Datname = $InputFile ;
    $Datname =~ s/\.txt$/.dat/ ;
    my $MaxPeak = 80 ;
    open( DAT , '>' , $Datname ) or print_error("Cannot open $Datname") ;
    for ( my $OneY = 0 ; $OneY < $MaxPeak ; $OneY++ ) {
	my $OutY = at($Matrix, $OneY, $OneX) ;
	$OutY = abs($OutY) ;
	my $OutX = abs(($MaxPeak - $OneY)* $EM::GridWidth * 0.001) ; 
	print DAT "$OutX $OutY\n" ;
    }
    close DAT ;

    ##################################################################



    $Mittelwert = $Mittelwert / ($EM::MatrixHeight * $EM::MatrixWidth) ;


    print "Average over all values: $Mittelwert \n";


	    print "Max: $Max \n" ;
	    print "Min: $Min \n" ;


    my $Check = calc_biot_savart(300) ;
    print "Biot-Savart for 300 nm = $Check [N/Am]\n" ;




    for ( my $Y = 0 ; $Y < $EM::MatrixHeight ; $Y++ ) {

	for ( my $X = 0 ; $X < $EM::MatrixWidth ; $X++ ) {
	    
#	    print "X = $X \nY = $Y \n" ;

#############################################################
# How to normalize the Matrix!?!

	    # Reset yourself to 256 Colors:

	    # my $NewVal = at($Matrix, $Y, $X) + (-$Min);
	    # $NewVal = $NewVal * 255 / (-$Min + $Max) ;



	    # Don't care for + or - Bz

	    my $NewVal = abs(at($Matrix, $Y, $X)) ;


	    # Set the average to 128 (+/- 2/4 * average)

#	    $NewVal = (2 * $Mittelwert) if abs(at($Matrix, $Y, $X)) > (2 * $Mittelwert) ;

#	    $NewVal = (4 * $Mittelwert) if abs(at($Matrix, $Y, $X)) > (4 * $Mittelwert) ;


	    # Set all values above A to A

#	    $NewVal = $EM::BSThickness if abs(at($Matrix, $Y, $X)) > $EM::BSThickness ;


	    # Set all changes to Matrix:

#	    $NewVal = $NewVal + $Mittelwert ;

	    set ($Matrix, $Y, $X, $NewVal);

	}
    }





    # 1st: show the Matrix on Screen (dev /XSERVE)

    dev '/XSERVE', {Recording => 1} ;

    ctab( lut_data("heat") ) ;
    imag($Matrix);
    hold ;
    cont($Matrix);


    # 2nd: redo all to a gif

    my $state = retrieve_state() ;
    dev '/GIF' ;
    replay $state ; 


    # 3rd: rename the gif

    my $BaseName = $InputFile ;
    $BaseName =~ s/\..*$// ;
    rename ("pgplot.gif" , "${BaseName}.gif") ;



# Output the Conducting lines into a gif

    if ( @AllLines ) {

	# create new Matrix:

	my $LineMatrix = zeroes( $EM::MatrixHeight, $EM::MatrixWidth );	

	# Print the conducting lines (stored in @AllLines)
	
	while ( @AllLines ) {
	    my $Lx = shift @AllLines ;
	    my $Ly = shift @AllLines ;
	    my $LineValue = 1 ;
	    
	    set ($LineMatrix, $Ly, $Lx, $LineValue);
	}

	dev '/GIF' ;
	imag($LineMatrix) ;
	rename ("pgplot.gif" , "${BaseName}line.gif") ;

    }






}