#!/usr/bin/perl

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

use strict;
use Sci qw($pi2 $e0 $pi $kB $e $e0 $me $h $hbar $KToeV erfc);
use Sci::OpticalMaterial;
use Sci::Material;
use Sci::Optimize;
use IniFile;

my $Header = 'ZnO-SCAM';
$Header = $ARGV[0] if(defined $ARGV[0]);

my $DataPath  = $Header . ".csv";
my $PrmPath   = $Header . ".prm";
my $OutputCSV = $Header . "-Simplex.csv";
print "Data: $DataPath\n";
print "Prm : $PrmPath\n";
print "Out : $OutputCSV\n";

my $ini = new IniFile(undef, 0, 0);
&ReadIni($PrmPath, $ini);

my $csv = new CSV();
$csv->Read($DataPath, 1) or die "$!\n";
my $pT     = $csv->GetXData(".*T.*");
my $pmu    = $csv->GetXData(".*mu.*");
my $pNe    = $csv->GetYData(".*Ne.*");
my $pSigma = $csv->GetYData(".*sigma.*");
my $nData  = @$pT;
print "nData=$nData\n";
#print "p=$pT, $pmu, $pNe, $pSigma\n";

my $out = new JFile($OutputCSV, "w") or die "$!: can not write to [$OutputCSV]";
my $material = new Material;

my $UseNeApproximation = (defined $ini->{UseNeApproximation})? $ini->{UseNeApproximation} : 0;
my $wNe = (defined $ini->{wNe})? $ini->{wNe} : 600;

#For Hall 
my $Opt = new Optimize; # For Hall effect
my $Method = "Amoeba::Simplex";
my $EPS = 1.0e-7;
my $nMaxIter = (defined $ini->{nMaxIter})? $ini->{nMaxIter} : 1000;

# EF計算収束条件
my $EFEPS    = 1.0e-5;
my $EFSEPS   = 1.0e-12;
#my $nMaxIter = 1000;
my $IsPrint  = 0;

# バンド端構造
my $Effme = 0.27;
my $Effmh = 1.0;
my $Ec    = 3.37;
my $Ev    = 0.0;
my $Eg    = $Ec - $Ev;
my $CBMultiplicity = 1.0;
my $VBMultiplicity = 1.0;

#my @Guess = (   $ND,  $ED, $d2,   $Ne2, $mu2, $Egb, $p[0], $u0[0], $p[1], $u0[1], $p[2], $u0[2]);
#S2 = 1.78986
# ドーパント
my $ND = (defined $ini->{ND})? $ini->{ND} : 6.545e14; #1.0e17;
my $ED = (defined $ini->{ED})? $ini->{ED} : 3.27886;  #$Ec - 0.030; # eV
my $NA = 0.0e18;
my $EA = $Ev + 0.10; # eV

# 層構造定義
my $FilmThickness = 300.0; # nm
my $d2  = (defined $ini->{d2})?  $ini->{d2}  : 0.3384;   #10.0;  # nm
my $Ne2 = (defined $ini->{Ne2})? $ini->{Ne2} : 1.705e18; #5e17;  # cm-3
my $mu2 = (defined $ini->{mu2})? $ini->{mu2} : 10.02;    #8.0;   # cm2/Vs

# バルク移動度定義
my @p  = ();
my @u0 = ();
my $Egb = (defined $ini->{Egb})? $ini->{Egb} : 0.0256; #0.032; # eV
$p[0]  = (defined $ini->{p0})?  $ini->{p0}  : -1.5;
$u0[0] = (defined $ini->{u00})? $ini->{u00} : 1826108;
$p[1]  = (defined $ini->{p1})?  $ini->{p1}  : 0;
$u0[1] = (defined $ini->{u01})? $ini->{u01} : 2789057;
$p[2]  = (defined $ini->{p2})?  $ini->{p2}  :  1.5;
$u0[2] = (defined $ini->{u02})? $ini->{u02} : 2163952;
$p[3]  = (defined $ini->{p3})?  $ini->{p3}  : -1.0;
$u0[3] = (defined $ini->{u03})? $ini->{u03} : 2163952;
$p[4]  = (defined $ini->{p4})?  $ini->{p4}  :  1.0;
$u0[4] = (defined $ini->{u04})? $ini->{u04} : 2163952;

# 伝導帯裾状態
#my $Eu    = 0.01;   # eV
#my $Ntail = 0.0e19; # cm-3
#my $Nt0   = 0.0;    # $Ntail / $Eu;
#my $Eth   = $Ec + 1;
my $Ecmin = $Ec - 0.0; # Minimum energy for CB integration
#printf "Ntail=%10.4g Nt0=%10.4g cm-3 E0=%6.3f eV\n", $Ntail, $Nt0, $Eu;

print("Initial values\n");
printf("UseNeApproximation=$UseNeApproximation\n");
printf("wNe=$wNe\n");
printf("  nMaxIter =$nMaxIter\n");
printf("  S2 =$ini->{S2}\n");
printf("  ND =%e cm-3\n", $ND);
printf("  ED =$ED eV\n");
printf("  d2 =$d2 nm\n");
printf("  Ne2=%e cm-3\n", $Ne2);
printf("  mu2=$mu2 cm2/Vs\n");
printf("  Egb=$Egb eV\n");
for(my $i = 0 ; $i < @p ; $i++) {
	print("  p$i=$p[$i]   u0$i=$u0[$i]\n");
}

my $T      = 300.0;

my $Nc   = $material->Nc($Effme*$me,  $CBMultiplicity, $T, 0) * 1.0e-6;
my $Dc0  = $material->Dc0($Effme*$me, $CBMultiplicity, 0) * $e**1.5 * 1.0e-6; # cm^-3/eV^1.5
my $Nv   = $material->Nv($Effmh*$me,  $VBMultiplicity, $T, 0) * 1.0e-6;
my $Dv0  = $material->Dv0($Effmh*$me, $VBMultiplicity, 0) * $e**1.5 * 1.0e-6; # cm^-3/eV^1.5
my $Ni   = $material->Ni($T, $Eg);
my $NeBase = 0.0;

print  "\n";
print  "Ec=$Ec eV\n";
print  "Ev=$Ev eV\n";
print  "ND=$ND cm-3  ED=$ED eV\n";
printf "Nc=%10.4g cm-3  Dc0=%10.4g cm-3/eV1.5\n", $Nc, $Dc0;
printf "Nv=%10.4g cm-3  Dv0=%10.4g cm-3/eV1.5\n", $Nv, $Dv0;
printf "Ni=%10.4g cm-3\n", $Ni;

my $dE     = $T * $KToeV * 0.1; #0.02;
my $Erange = $T * $KToeV * 7.0 + $Ec - $Ecmin;
my $nE     = int($Erange / $dE + 1.001);
$nE = int($nE / 2) * 2 + 1;
printf "Integration range: %6.3f, %10.3g eV step [%d]\n", $Erange, $dE, $nE;

my (@Energyc, @Dc, @Energyv, @Dv);
&CalDOSc();
&CalDOSv();

my ($EFopt, $S) = &CalEF($T, $NeBase, $EFEPS, $EFSEPS, 0);
my $Ne  = &CalNe($EFopt, $T);
my $Nh  = &CalNh($EFopt, $T);
printf "EF=%g eV: Ne=%7.4g Nh=%7.4g  (S=%5.3g)\n", $EFopt, $Ne, $Nh, $S;


for(my $i = 0 ; $i < $nData ; $i++) {
	my $T = $pT->[$i];
	my ($u, $n) = &CalMobility($T, $FilmThickness, $d2, $ND, $ED, $d2, $Ne2, $mu2, $Egb, \@p, \@u0);
	printf "%d: T=%4.3f  ucal=%8.3g (%8.3g)  Necal=%8.3g (%8.3g)\n",
			$i, $T, $u, $pmu->[$i], $n, $pNe->[$i];
}

my @Guess = (   $ND,  $ED, $d2,   $Ne2, $mu2, $Egb, $p[0], $u0[0], $p[1], $u0[1], $p[2], $u0[2], $p[3], $u0[3], $p[4], $u0[4]);
#my @OptId = (     0,    0,   0,      0,    0,    0,     0,      1,     0,      1,     0,      1,     0,      1,     0,      1);
#my @OptId = (     0,    0,   0,      0,    0,    0,     0,      1,     0,      1,     0,      1,     0,      1,     0,      1);
my @OptId = (     1,    1,   0,      1,    1,    0,     0,      1,     0,      1,     0,      1,     0,      1,     0,      1);
my @Scale = (3.0e17, 0.03,  10, 1.0e20,    8, 0.01,   0.1,  1.0e8,   0.1,  1.0e7,   0.1,  1.0e2,   0.1,  1.0e5,   0.1,  1.0e6);
my $iPrintLevel = 0;

my ($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				\&MinimizationFunc, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";

($ND, $ED, $d2, $Ne2, $mu2, $Egb, $p[0], $u0[0], $p[1], $u0[1], $p[2], $u0[2], $p[3], $u0[3], $p[4], $u0[4]) = @$OptVars;

$out->print("T(K),mu(obs),mu(cal),Ne(obs),Ne(cal),mu(bulk),Ne(bulk)\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $T = $pT->[$i];
	my ($u, $n, $u1, $n1) = &CalMobility($T, $FilmThickness, $d2, $ND, $ED, $d2, $Ne2, $mu2, $Egb, \@p, \@u0);

	printf "%d: T=%4.3f  ucal=%8.3g (%8.3g)  Necal=%8.3g (%8.3g)\n",
			$i, $T, $u, $pmu->[$i], $n, $pNe->[$i];
	$out->print("$T,$pmu->[$i],$u,$pNe->[$i],$n,$u1,$n1\n");
}
$out->print("S2,$MinVal,\n");
$out->print("ND,$ND,cm-3\n");
$out->print("ED,$ED,eV\n");
$out->print("d2,$d2,nm\n");
$out->print("Ne2,$Ne2,cm-3\n");
$out->print("mu2,$mu2,cm2/Vs\n");
$out->print("Egb,$Egb,eV\n");
for(my $i = 0 ; $i < @p ; $i++) {
	$out->print("p$i,$p[$i],,u0$i,$u0[$i]\n");
}

$out->print("T(K),mu(cal),Ne(cal),mu(bulk),Ne(bulk)\n");
for(my $T = 30 ; $T < 350 ; $T += 10) {
	my ($u, $n, $u1, $n1) = &CalMobility($T, $FilmThickness, $d2, $ND, $ED, $d2, $Ne2, $mu2, $Egb, \@p, \@u0);

	printf "  T=%4.3f  ucal=%8.3g  Necal=%8.3g\n", $T, $u,$n;
	$out->print("$T,$u,$n,$u1,$n1\n");
}
$out->Close();

print("S2=$MinVal\n");
printf("  ND=%e cm-3\n", $ND);
printf("  ED=$ED eV\n");
printf("  d2=$d2 nm\n");
printf("  Ne2=%e cm-3\n", $Ne2);
printf("  mu2=$mu2 cm2/Vs\n");
printf("  Egb=$Egb eV\n");
for(my $i = 0 ; $i < @p ; $i++) {
	print("  p$i=$p[$i]   u0$i=$u0[$i]\n");
}

open(OUT, ">$PrmPath") or die "$!\n";
printf(OUT "UseNeApproximation=$UseNeApproximation\n");
printf(OUT "wNe=$wNe\n");
printf(OUT "nMaxIter=$nMaxIter\n");
printf(OUT "S2=$MinVal\n");
printf(OUT "ND=%e\n", $ND);
printf(OUT "ED=$ED\n");
printf(OUT "d2=$d2\n");
printf(OUT "Ne2=%e\n", $Ne2);
printf(OUT "mu2=$mu2\n");
printf(OUT "Egb=$Egb\n");
for(my $i = 0 ; $i < @p ; $i++) {
	printf(OUT "p$i=$p[$i]\n");
	printf(OUT "u0$i=$u0[$i]\n");
}
close(OUT);

exit;

sub ReadIni
{
	my ($PrmPath, $ini) = @_;
	open(IN, $PrmPath) or return;
	while(1) {
		my $line = <IN>;
		last if(!defined $line);
		my ($key, $val) = ($line =~ /^\s*(.*?)\s*=\s*(.*)\s*$/);
		next if(!defined $val);
#print "$key: $val\n";
		$ini->{$key} = $val;
	}
	close(IN);
	return $ini;
}

sub CalMobility
{
	my ($T, $FilmThickness, $d2, $Nd, $Ed, $d2, $Ne2, $mu2, $Egb, $pp, $pu0) = @_;

	my $d1 = $FilmThickness  - $d2;
	$ND = $Nd;
	$ED = $Ed;

	my $ET = $kB * $T / $e; # eV;
#print "$i: T=$T   E=$ET eV\n";

if($UseNeApproximation) {
	$Nc     = $material->Nc($Effme*$me,  $CBMultiplicity, $T, 0) * 1.0e-6;
#	$Nv     = $material->Nv($Effmh*$me,  $VBMultiplicity, $T, 0) * 1.0e-6;
}
else {
#	$Nc     = $material->Nc($Effme*$me,  $CBMultiplicity, $T, 0) * 1.0e-6;
#	$Nv     = $material->Nv($Effmh*$me,  $VBMultiplicity, $T, 0) * 1.0e-6;
	$dE     = $T * $KToeV * 0.05; #0.1 #0.02;
	$Erange = $T * $KToeV * 10.0 + $Ec - $Ecmin;
	$nE     = int($Erange / $dE + 1.001);
	$nE = int($nE / 2) * 2 + 1;
	&CalDOSc();
	&CalDOSv();
}

	my $invu = 0.0;
	for(my $j = 0 ; $j < @$pp ; $j++) {
		if($pp->[$j] == 0.0) {
			$invu += 1.0 / $pu0->[$j];
		}
		else {
			$invu += 1.0 / ($pu0->[$j] * $T ** $pp->[$j]);
		}
	}
	my $mu1 = exp(-$Egb / $ET) / $invu;
	my $Ne1;
	if($UseNeApproximation) {
		$Ne1 = sqrt(abs($Nc*$ND)) * exp(-($Ec-$ED) / $ET / 2.0);
	}
	else {
		my ($EFopt, $S) = &CalEF($T, $NeBase, $EFEPS, $EFSEPS, 0);
		$Ne1  = &CalNe($EFopt, $T);
	}
	my $u = ($mu1*$mu1*$Ne1*$d1 + $mu2*$mu2*$Ne2*$d2)    / ($mu1     *$Ne1*$d1 + $mu2     *$Ne2*$d2);
	my $n = ($mu1     *$Ne1*$d1 + $mu2     *$Ne2*$d2)**2 / ($mu1*$mu1*$Ne1*$d1 + $mu2*$mu2*$Ne2*$d2) / $FilmThickness;
	return ($u, $n, $mu1, $Ne1);
}

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub MinimizationFunc {
	my ($pVars, $iPrintLevel) = @_;
	my ($ND, $ED, $d2, $Ne2, $mu2, $Egb, @p, @u0);
	($ND, $ED, $d2, $Ne2, $mu2, $Egb, $p[0], $u0[0], $p[1], $u0[1], $p[2], $u0[2], $p[3], $u0[3], $p[4], $u0[4]) = @$pVars;

	my $S = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $T = $pT->[$i];
		my ($u, $n) = &CalMobility($T, $FilmThickness, $d2, $ND, $ED, $d2, $Ne2, $mu2, $Egb, \@p, \@u0);
#print "$i: ucal=$u  uobs=$pmu->[$i]\n";
		my $du = $u - $pmu->[$i];
		my $dNe = ($n > 0.0)? log($n) - log($pNe->[$i]) : -log($pNe->[$i]);
		$S += $du*$du + $wNe * $dNe*$dNe;
		$S += 1.0e10 if($u0[0] < 0.0);
		$S += 1.0e10 if($u0[1] < 0.0);
		$S += 1.0e10 if($u0[2] < 0.0);
		$S += 1.0e10 if($u0[3] < 0.0);
		$S += 1.0e10 if($u0[4] < 0.0);
		$S += 1.0e10 if($d2    < 0.0 or $d2 > $FilmThickness);
		$S += 1.0e10 if($ND    < 0.0);
		$S += 1.0e10 if($mu2   < 0.0);
##		$S += 1.0e10 if($Ne2   < 0.0);
##		$S += 1.0e10 if($n     < 0.0);
	}
	$S = sqrt($S/2.0/$nData);
printf "Vars: %8.3g, %5.3g, %5.3g, %8.3g, %5.3g, %8.3g, %5.2g, %8.3g, %5.2g, %8.3g, %5.2g, %8.3g, %5.2g, %8.3g, %5.2g, %8.3g: S=%8.4f\n", 
		$ND, $ED, $d2, $Ne2, $mu2, $Egb, $p[0], $u0[0], $p[1], $u0[1], $p[2], $u0[2], $p[3], $u0[3], $p[4], $u0[4], $S;
	return $S;
}

sub PDLDisplayFunc {
	my ($piddle) = @_;
#print "Display: piddle=$piddle\n";
}

#===================================
# Material Module
#===================================
# 裾状態を含んだパラボリック伝導帯の状態密度
sub CalDOSc
{
#	my (@Dtail);

#	my $IsTail = 1;
#	$IsTail = 0 if($Nt0 == 0.0);
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Ecmin + $i * $dE;
		$Energyc[$i] = $E;
		$Dc[$i] = ($E > $Ec)? $Dc0 * sqrt($E - $Ec) : 0.0;
#		$Dtail[$i] = $Nt0 * exp(($E-$Eth) / $Eu);
#		if($IsTail and $Dtail[$i] > $Dc[$i]) {
#			$Dc[$i] = $Dtail[$i];
#		}
#		elsif($IsTail) {
#			my $Ntrap = $Nt0 * $Eu * exp(($E-$Eth) / $Eu);
#print "Eborder=$E  Eth=$Eth eV  Ntrap=$Ntrap Nt0=$Nt0 cm-3\n";
#			$IsTail = 0;
#		}
	}

#	my $f = "Dc.csv";
#	my $out = new JFile($f, "w");
#	$out->print("E,Dc,Dtail\n");
#	for(my $i = 0 ; $i < $nE ; $i++) {
#		$out->print("$Energyc[$i],$Dc[$i],$Dtail[$i]\n");
#	}
#	$out->Close();
}

# パラボリック価電子帯の状態密度
sub CalDOSv
{
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Ev - $i * $dE;
		$Energyv[$i] = $E;
		$Dv[$i] = ($E < $Ev)? $Dv0 * sqrt($Ev - $E) : 0.0;
	}
}

# 二分法でフェルミ準位を計算
sub CalEF
{
	my ($T, $NeBase, $xEPS, $SEPS, $IsPrint) = @_;

my $nMaxIter=1000;
	my $opt = new Optimize; # For EF
	my ($pVars, $S) = $opt->BinarySearch(
							$Ev-1.0, $Ec+1.0, sub { &CalS($T, $NeBase, @_); },
							$nMaxIter, $xEPS, $SEPS, $IsPrint
						);
	return ($pVars->[0], $S);
}

# 過剰電荷の計算。過剰電荷＝＝0のEFを探すために使う
sub CalS 
{
	my ($T, $NeBase, $EF, $iPrintLevel) = @_;

	my $Ne  = &CalNe($EF, $T);
	my $Nh  = &CalNh($EF, $T);
	my $NDp = &CalNDp($EF, $T);
	my $NAm = &CalNAm($EF, $T);
#print "EF=$EF  Ne=$Ne ND=$ND (NDp=$NDp)  Nh=$Nh  NA=$NA(NAm=$NAm)\n";
	my $dN = $Ne - $NDp - ($Nh - $NAm) - $NeBase;
	my $S = $dN;
printf("EF=%6.3f: Ne=$Ne Nh=$Nh (%g)\n", $EF, $S) if($iPrintLevel >= 2);
	return $S;
}

# 伝導帯中の占有電子数
# EDが伝導帯中にある場合でも、ドナー準位を占める電子数は数えていない
sub CalNe
{
	my ($EF, $T) = @_;
	my @Dfc;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Energyc[$i];
		$Dfc[$i] = $Dc[$i] * &FEe($E, $EF, $T);
	}
	return  Algorism::SinglePointIntegrateBySimpson(\@Energyc, \@Dfc);
}

# 価電子中の占有正孔数
# EAが価電子帯中にある場合でも、アクセプター準位を占める正孔数は数えていない
sub CalNh
{
	my ($EF, $T) = @_;
	my @Dfv;
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Energyv[$i];
		$Dfv[$i] = $Dv[$i] * &FEh($E, $EF, $T);
	}
	return -Algorism::SinglePointIntegrateBySimpson(\@Energyv, \@Dfv);
}

# イオン化ドナー密度
sub CalNDp
{
	my ($EF, $T) = @_;
	return $ND * &FEh($ED, $EF, $T);
}

# イオン化アクセプター密度
sub CalNAm
{
	my ($EF, $T) = @_;
	return $NA * &FEe($EA, $EF, $T);
}

# 電子のFermi-Dirac分布のエネルギー微分
sub dFEedE
{
	my ($E, $EF, $T) = @_;
	return 0.0 if($T == 0.0);

	my $ex  = exp(($E-$EF) / $T / $KToeV);
	my $ex1 = 1.0 + $ex;
	return -1.0 / $T / $KToeV * $ex / ($ex1*$ex1);
}

# 電子のFermi-Dirac分布
sub FEe
{
	my ($E, $EF, $T) = @_;
	if($T == 0.0) {
		return 0.0 if($E > $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($E-$EF) / $T / $KToeV) );
}

# 正孔のFermi-Dirac分布
sub FEh
{
	my ($E, $EF, $T) = @_;
	if($T == 0.0) {
		return 0.0 if($E < $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($EF-$E) / $T / $KToeV) );
}