#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
print "\n\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC);
}

use strict;
#use warnings;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci::Science;

my $e0  = Sci::e0();
my $e   = Sci::e();
my $pi  = Sci::pi();
my $pi2 = 2.0 * $pi;
my $c   = Sci::c();

my $OutFile = "fs.csv";
my $OutFile2 = "fs2.csv";

my $lambda  = 800e-9;  # m 波長
my $dtPulse = 10e-15; # s  パルス幅
my $LPulse  = $dtPulse * $c; # m パルス長
my $k0      = $pi2 / $lambda;

#for charp
#my $deltaPhi = 5.0 * $pi/4.0 / 30.0; #移送ずれ係数 rad/nm, Δφ=$dPhi * (Lambda / $lambda - 1.0);
my $deltaPhi = 0.0;

my $a  = $LPulse; # m
my $x0 = -5.0 * $LPulse;
my $x1 = 10.0 * $LPulse;
my $nx = 3001;

my $kstart = -30.0 / $a;
my $kend   =  30.0 / $a;
my $nk     = 1001;

#my $kFTmin = 1.0 * $kstart;
#my $kFTmax = 1.0 * $kend;
my $kFTmin = 0.02 * $kstart;
my $kFTmax = 0.02 * $kend;
my $nkFT   = 11;
my $MaxLambda = 0.0; #100; # um

my $xstep   = ($x1 - $x0) / ($nx - 1);
my $kstep   = ($kend - $kstart) / ($nk - 1);
my $kFTstep = ($kFTmax - $kFTmin) / ($nkFT - 1);

print "x(OUT): $x0 - $x1, $xstep ($nx)\n";
print "k(OUT): $kstart - $kend, $kstep ($nk)\n";
print "k(FT) : $kFTmin - $kFTmax, $kFTstep ($nkFT)\n";

open(OUT, ">$OutFile") or die "$!\n";
print "k(cm-1),lambda(nm),A(cos kx),B(sin kx)\n";
print OUT "k(cm-1),lamda(nm),A(cos kx),B(sin kx),|R|\n";

for(my $i = 0 ; $i < $nk ; $i++) {
	my $k = $kstart + $i * $kstep + $k0;
	my $kp = $k + $k0;
	my $km = $k - $k0;
	my $l;
	if($k == 0.0) {
		$l = 1e9;
	}
	else {
		$l = $pi2 / $k * 1.0e9;
	}
	my ($A, $B);
	if($km == 0.0) {
		$A = (sin($kp*$a)/$kp + $a) / 2.0 / $pi2;
		$B = ((1.0 - cos($kp*$a))/$kp) / 2.0 / $pi2;
	}
	else {
		$A = (sin($kp*$a)/$kp + sin($km*$a)/$km) / 2.0 / $pi2;
		$B = ((1.0 - cos($kp*$a))/$kp + (1.0 - cos($km*$a))/$km) / 2.0 / $pi2;
	}
printf("%8.4g,%12.6g,%12.6g\n", $l, $A, $B);
printf(OUT "%8.4g,%8.4g,%12.6g,%12.6g,%12.6g\n", $k*100, $l, $A, $B, sqrt($A*$A+$B*$B));
}

open(OUT, ">$OutFile2") or die "$!\n";
print "x(m),f\n";
print OUT "x(um),f";
for(my $j = 0 ; $j < $nkFT ; $j++) {
	my $k = $kFTmin + $j * $kFTstep + $k0;
	my $l;
	if($k == 0.0) {
		$l = 1e9;
	}
	else {
		$l = $pi2 / $k * 1.0e9;
	}
	next if($MaxLambda > 0.0 and abs($l) > $MaxLambda);
	printf OUT ",A(%6.3gnm)", $l;
}
print OUT "\n";

for(my $i = 0 ; $i < $nx ; $i++) {
	my $x = $x0 + $i * $xstep;
	my $f = 0.0;
	my (@k, @A, @B);
	for(my $j = 0 ; $j < $nkFT ; $j++) {
		my $k = $kFTmin + $j * $kFTstep + $k0;
		my $kp = $k + $k0;
		my $km = $k - $k0;
		$k[$j] = $k;
		my $l;
		if($k == 0.0) {
			$l = 1e6;
		}
		else {
			$l = $pi2 / $k * 1.0e6;
		}
		next if($MaxLambda > 0.0 and abs($l) > $MaxLambda);
		if($km == 0.0) {
			$A[$j] = (sin($kp*$a)/$kp + $a) / 2.0 / $pi2;
			$B[$j] = ((1.0 - cos($kp*$a))/$kp) / 2.0 / $pi2;
		}
		else {
			$A[$j] = (sin($kp*$a)/$kp + sin($km*$a)/$km) / 2.0 / $pi2;
			$B[$j] = ((1.0 - cos($kp*$a))/$kp + (1.0 - cos($km*$a))/$km) / 2.0 / $pi2;
		}
		my $dPhi = $deltaPhi * ($l*1e-3 / $lambda - 1.0);
#print "dp=$dPhi (", $l/$lambda*1e-3, ")[$l]\n";
		$f += $A[$j] * cos($k*$x + $dPhi) + $B[$j] * sin($k*$x + $dPhi);
	}
	$f *= $kFTstep;
printf("%8.4g,%12.6g\n", $x*1.0e6, $f);
printf(OUT "%8.4g,%12.6g", $x*1.0e6, $f);
	for(my $j = 0 ; $j < $nkFT ; $j++) {
		my $l;
		if($k[$j] == 0.0) {
			$l = 1e6;
		}
		else {
			$l = $pi2 / $k[$j] * 1.0e6;
		}
		next if($MaxLambda > 0.0 and abs($l) > $MaxLambda);
		printf OUT ",%12.6g", $A[$j]*cos($k[$j]*$x)*$kFTstep + $B[$j]*sin($k[$j]*$x)*$kFTstep;
	}
	print OUT "\n";
}

close(OUT);