Start implementing rawfield and tests

wlmb · wlmb · commit fbfac241fd61 · 2026-01-02T16:21:47.000-06:00
diff --git a/Changes b/Changes
@@ -2,6 +2,9 @@ Revision history for Photonic
 
 {{$NEXT}}
 
+  - Program rawfield to obtain the unnormalized fields, assuming the
+    external source is an external polarization. Currently, modified
+    LE::NR2::Field, LE::S::Field, Roles::Field, WE::S::Field, and the tests.
   - Add conversion routines between frequency, energy and wavelength
   - Avoid duplicate arguments in cgtsv
 
diff --git a/lib/Photonic/LE/NR2/Field.pm b/lib/Photonic/LE/NR2/Field.pm
@@ -86,6 +86,7 @@ evaluation of the field
 use namespace::autoclean;
 use PDL::Lite;
 use PDL::NiceSlice;
+use PDL::Constants qw(PI);
 use Photonic::LE::NR2::Haydock;
 use Photonic::Utils qw(cgtsv GtoR linearCombineIt);
 use Photonic::Types -all;
@@ -165,15 +166,16 @@ sub _build_field {
 
 sub _build_rawfield {
     my $self=shift;
-    my $Es = $self->_Fn*$self->u/$self->epsA;
+    my $Es = -4*PI*$self->_Fn*$self->u/$self->epsA; # assume drive is external polarization
     my $states=$self->haydock->states->dummy(0);
     my $nrGnorm = $self->haydock->GNorm;
     #field is cartesian,nx,ny...
     my $field_G=linearCombineIt($Es, $nrGnorm*$states); #En ^G|psi_n>
     $field_G *= $self->filter->(*1) if $self->has_filter;
     # fourier transform vector field.
     my $field_R=GtoR($field_G, $self->haydock->ndims, 1);
-    $field_R*=$self->haydock->B->nelem; #scale FFT
+    my $b0=$self->haydock->bs->((0));
+    $field_R*=$b0*$self->haydock->B->nelem; #scale FFT
     return $field_R; #result is cartesian, nx, ny,...
 }
 
diff --git a/lib/Photonic/LE/S/Field.pm b/lib/Photonic/LE/S/Field.pm
@@ -72,6 +72,7 @@ evaluation of the field
 use namespace::autoclean;
 use PDL::Lite;
 use PDL::NiceSlice;
+use PDL::Constants qw(PI);
 use Photonic::LE::S::Haydock;
 use Photonic::Utils qw(cgtsv GtoR linearCombineIt);
 use Photonic::Types -all;
@@ -95,8 +96,10 @@ sub _build_epsL {
     $self->epsL; # danger of infinite recursion?
 }
 
+has '_Fn'=>(is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+	    documentation=>'Solution of (epsLL)F=(1,0...) in Haydock basis');
 
-sub _build_field {
+sub _build__Fn{
     my $self=shift;
     my $as=$self->haydock->as;
     my $bs=$self->haydock->bs;
@@ -118,15 +121,22 @@ sub _build_field {
     # Add spinor normalization.
     my $epsL=1/$result->((0));
     $self->_epsL($epsL);
+    return $result;
+}
+
+sub _build_field {
+    my $self=shift;
+    my $result=$self->_Fn;
+    my $epsL=$self->epsL;
     my $norm=sqrt(2)*$epsL;
     # Normalize result so macroscopic field is 1.
     my $Es = $result*$norm;
     #states are xy,nx,ny...
-    my $stateit=$self->haydock->states->slice("*1");
+    my $states=$self->haydock->states->slice("*1"); #dummy(0)
     #pmGnorm is xy,pm,nx,ny...
     my $pmGNorm=$self->haydock->pmGNorm;
     #field is xy,pm,nx,ny...
-    my $field_G=linearCombineIt($Es, $pmGNorm*$stateit); #En ^G|psi_n>
+    my $field_G=linearCombineIt($Es, $pmGNorm*$states); #En ^G|psi_n>
     #Choose +k
     my $Esp=$field_G->(:,(0)); #xy,nx,ny
     $Esp *= $self->filter->(*1) if $self->has_filter;
@@ -136,6 +146,28 @@ sub _build_field {
     return $field_R; #result is xy,nx,ny,...
 }
 
+sub _build_rawfield {
+    my $self=shift;
+    my $Es=-4*PI*$self->_Fn; # Drive is external polarization
+    # Multiply by sqrt(2)? Size of firststate is in b0? Included in $Es?
+    #states are xy,nx,ny...
+    my $states=$self->haydock->states->slice("*1"); #dummy(0)
+    #pmGnorm is xy,pm,nx,ny...
+    my $pmGNorm=$self->haydock->pmGNorm;
+    #field is xy,pm,nx,ny...
+    my $field_G=linearCombineIt($Es, $pmGNorm*$states); #En ^G|psi_n>
+    #Choose +k
+    my $Esp=$field_G->(:,(0)); #xy,nx,ny
+    $Esp *= $self->filter->(*1) if $self->has_filter;
+    #get cartesian out of the way, fourier transform, put cartesian.
+    my $field_R=GtoR($Esp, $self->haydock->B->ndims, 1);
+    $field_R*=$self->haydock->B->nelem; #scale FFT
+    my $b0=$self->haydock->bs->slice((0));    # Should I multiply by $b0?
+    $field_R*=$b0*sqrt(2); # sqrt 2 from spinor normalization
+    return $field_R; #result is xy,nx,ny,...
+}
+
+
 __PACKAGE__->meta->make_immutable;
 
 1;
diff --git a/lib/Photonic/Roles/Field.pm b/lib/Photonic/Roles/Field.pm
@@ -96,6 +96,7 @@ real space field normalized to macroscopic field in format cartesian, nx, ny,...
 =item * rawfield
 
 real space field normalized to external source in format cartesian, nx, ny,... (lazy-built).
+Source is interpreted as external polarization.
 
 =item * nh
 
diff --git a/lib/Photonic/WE/S/Field.pm b/lib/Photonic/WE/S/Field.pm
@@ -61,15 +61,16 @@ Consumes L<Photonic::Roles::Field>
 use namespace::autoclean;
 use PDL::Lite;
 use PDL::NiceSlice;
+use PDL::Constants qw(PI);
 use Photonic::WE::S::Haydock;
 use Photonic::Utils qw(cgtsv GtoR linearCombineIt);
 use Photonic::Types -all;
 use Moo;
 use MooX::StrictConstructor;
 
 # Temporary:
-has 'rawfield'=>(is=>'lazy', isa=>PDLComplex,
-           documentation=>'Calculated real space field, unnormalized');
+#has 'rawfield'=>(is=>'lazy', isa=>PDLComplex,
+#           documentation=>'Calculated real space field, unnormalized');
 
 with 'Photonic::Roles::Field';
 
@@ -105,14 +106,14 @@ sub _build_rawfield {
     my $Es=$self->haydock->applyMetric($field_G);
     #Comment as unnormalized
     #$Es*=$bs->((0))/$self->haydock->metric->epsilon;
-    my $Esp=$Es(:,(0)); # choose +k spinor component.
+    my $Esp=sqrt(2)*$Es(:,(0)); # choose +k spinor component.
     $Esp *= $self->filter->(*1) if $self->has_filter;
     ##get cartesian out of the way, fourier transform, put cartesian.
     my $field_R=GtoR($Esp, $ndims, 1);
     $field_R*=$self->haydock->B->nelem; #scale??
     my $b0=$self->haydock->bs->slice('(0)'); # #First state normalization factor
     $field_R*=$b0; #scale??
-    $field_R/=$self->haydock->B->nelem; #normalize FT?
+    $field_R*= -4*PI; # Interpreting the first state as external polarization.
     return $field_R; #result is xy,nx,ny...
 }
 
diff --git a/t/field-lenr2.t b/t/field-lenr2.t
@@ -31,6 +31,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA  02110-1301 USA
 use strict;
 use warnings;
 use PDL;
+use PDL::Constants qw(PI);
 use Photonic::LE::NR2::Haydock;
 use Photonic::LE::NR2::Field;
 use Photonic::LE::NR2::SHP;
@@ -55,8 +56,8 @@ my $fla=1/$ea;
 my $flb=1/$eb;
 my $fproml=$fla*(1-$gl->f)+$flb*($gl->f);
 my $flex=1/$fproml;
-($fla, $flb)=map {$_/$fproml} ($fla, $flb);
-my $flx=($fla*(1-$B)+$flb*$B)->slice("*1");
+my ($flan, $flbn)=map {$_/$fproml} ($fla, $flb);
+my $flx=($flan*(1-$B)+$flbn*$B)->slice("*1");
 ok(Cagree($flv, $flx), "1D long field") or diag "got: $flv\nexpected: $flx";
 ok(Cagree($fle, $flex), "1D long response") or diag "got: $fle\nexpected: $flex";
 #View 2D from 1D superlattice.
@@ -71,6 +72,12 @@ ok(Cagree($ftv, $ftx), "1D trans field") or diag "got: $ftv\nexpected: $ftx";;
 my $fpromt=$ea*(1-$gt->f)+$eb*($gt->f);
 ok(Cagree($fte, $fpromt), "1D trans response") or diag "got: $fte\nexpected: $fpromt";
 
+# check raw fields
+# Longitudinal
+my $flv_raw=$flo->rawfield;
+my $flx_raw=-4*PI*($fla*(1-$B)+$flb*$B)->dummy(0);
+ok(Cagree($flv_raw, $flx_raw), "1D long raw field");
+
 
 my ($dA, $dB) = (0, 1); # vacuum, then anything as is normalised to dB
 my $nrshp=Photonic::LE::NR2::SHP->new(
diff --git a/t/field-les.t b/t/field-les.t
@@ -31,15 +31,19 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA  02110-1301 USA
 use strict;
 use warnings;
 use PDL;
+use PDL::Constants qw(PI);
 use Photonic::LE::S::Haydock;
 use Photonic::LE::S::Field;
-
-use Test::More tests => 4;
+use Photonic::Geometry::FromB;
+use Photonic::Geometry::FromEpsilon;
+use Test::More tests => 6;
 use lib 't/lib';
 use TestUtils;
 
-my $ea=1+2*i;
-my $eb=3+4*i;
+#my $ea=1+2*i;
+#my $eb=3+4*i;
+my $ea=1+0*i;
+my $eb=3+0*i;
 #Check field for simple 1D system
 my $B=zeroes(11)->xvals<5; #1D system
 my $epsilon=$ea*(1-$B)+$eb*$B;
@@ -48,6 +52,7 @@ my $haydock=Photonic::LE::S::Haydock->new(geometry=>$gl, nh=>10,
    keepStates=>1, epsilon=>$epsilon);
 my $flo=Photonic::LE::S::Field->new(haydock=>$haydock, nh=>10);
 my $flv=$flo->field;
+
 my $fle=$flo->epsL;
 my $fla=1/$ea;
 my $flb=1/$eb;
@@ -57,7 +62,6 @@ my $flex=1/$fproml;
 my $flx=($fla*(1-$B)+$flb*$B)->slice("*1");
 ok(Cagree($flv, $flx), "1D long field") or diag "got: $flv\nexpected: $flx";
 ok(Cagree($fle, $flex), "1D long response") or diag "got: $fle\nexpected: $flex";
-
 #View 2D from 1D superlattice.
 my $Bt=zeroes(1,11)->yvals<5; #2D flat system
 my $epsilont=$ea*(1-$Bt)+$eb*$Bt;
@@ -71,3 +75,12 @@ my $ftx=pdl([1, 0])->r2C;
 ok(Cagree($ftv, $ftx), "1D trans field") or diag "got: $ftv\nexpected: $ftx";;
 my $fpromt=$ea*(1-$gt->f)+$eb*($gt->f);
 ok(Cagree($fte, $fpromt), "1D trans response") or diag "got: $fte\nexpected: $fpromt";
+
+# check raw fields
+my $flv_raw=$flo->rawfield;
+my $flx_raw=-4*PI*((1-$B)/$ea+$B/$eb)->dummy(0);
+ok(Cagree($flv_raw, $flx_raw), "1D long raw field");
+
+my $ftv_raw=$fto->rawfield;
+my $ftx_raw=-(4*PI+0*i)/($ea*(1-$gt->f)+$eb*$gt->f)*pdl([1,0]);
+ok(Cagree($ftv_raw, $ftx_raw), "1D transverse raw field");
diff --git a/t/field-wes.t b/t/field-wes.t
@@ -31,45 +31,92 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA  02110-1301 USA
 use strict;
 use warnings;
 use PDL;
+use PDL::Constants qw(PI);
 use Photonic::WE::S::Haydock;
 use Photonic::WE::S::Metric;
 use Photonic::WE::S::Field;
 
-use Test::More tests => 2;
+use Test::More tests => 4;
 use lib 't/lib';
 use TestUtils;
 
 my $ea=r2C(1);
 my $eb=3+4*i;
 
-#Check field for simple 1D system. Longitudinal case
-my $B=zeroes(11)->xvals<5; #1D system
-my $epsilon=$ea*(1-$B)+$eb*$B;
-my $gl=Photonic::Geometry::FromB->new(B=>$B); #long
-my $ml=Photonic::WE::S::Metric->new(geometry=>$gl, epsilon=>pdl(1),
-   wavenumber=>pdl(1), wavevector=>pdl([0.01]));
-my $haydock=Photonic::WE::S::Haydock->new(
-   metric=>$ml, nh=>10, keepStates=>1, polarization=>pdl([1])->r2C,
-   epsilon=>$epsilon);
-my $flo=Photonic::WE::S::Field->new(haydock=>$haydock, nh=>10);
-my $flv=$flo->field;
-my $fla=1/$ea;
-my $flb=1/$eb;
-my $fproml=$fla*(1-$gl->f)+$flb*($gl->f);
-($fla, $flb)=map {$_/$fproml} ($fla, $flb);
-my $flx=($fla*(1-$B)+$flb*$B)->transpose;
-ok(Cagree($flv, $flx), "1D long field");
+#Check field for simple 1D system.
+{
+    #Longitudinal case
+    my $B=zeroes(11)->xvals<5; #1D system
+    my $epsilon=$ea*(1-$B)+$eb*$B;
+    my $gl=Photonic::Geometry::FromB->new(B=>$B); #long
+    my $ml=Photonic::WE::S::Metric->new(geometry=>$gl, epsilon=>pdl(1),
+					wavenumber=>pdl(1), wavevector=>pdl([0.01]));
+    my $haydock=Photonic::WE::S::Haydock->new(
+	metric=>$ml, nh=>10, keepStates=>1, polarization=>pdl([1])->r2C,
+	epsilon=>$epsilon);
+    my $flo=Photonic::WE::S::Field->new(haydock=>$haydock, nh=>10);
+    my $flv=$flo->field;
+    my $fla=1/$ea;
+    my $flb=1/$eb;
+    my $fproml=$fla*(1-$gl->f)+$flb*($gl->f);
+    ($fla, $flb)=map {$_/$fproml} ($fla, $flb);
+    my $flx=($fla*(1-$B)+$flb*$B)->transpose;
+    ok(Cagree($flv, $flx), "1D long field");
+}
+
+{
+    #View 2D from 1D superlattice. Long wavelength transverse case
+    my $Bt=zeroes(1,11)->yvals<5; #2D flat system
+    my $epsilont=$ea*(1-$Bt)+$eb*$Bt;
+    my $gt=Photonic::Geometry::FromB->new(B=>$Bt); #trans
+    my $mt=Photonic::WE::S::Metric->new(geometry=>$gt, epsilon=>pdl(1),
+					wavenumber=>pdl(0.001), wavevector=>pdl([0,0.0001]));
+    my $nt=Photonic::WE::S::Haydock->new(
+	metric=>$mt, nh=>10, keepStates=>1, polarization=>pdl([1,0])->r2C,
+	epsilon=>$epsilont);
+    my $fto=Photonic::WE::S::Field->new(haydock=>$nt, nh=>10);
+    my $ftv=$fto->field;
+    my $ftx=r2C(pdl [1, 0]);
+    ok(Cagree($ftv, $ftx), "1D trans field");
+}
+
+#Check rawfields for simple 1D system.
+# Longitudinal case
+{
+    my $B=zeroes(11)->xvals<5; #1D system
+    my $epsilon=$ea*(1-$B)+$eb*$B;
+    my $gl=Photonic::Geometry::FromB->new(B=>$B); #long
+    my $ml=Photonic::WE::S::Metric->new(geometry=>$gl, epsilon=>pdl(1),
+					wavenumber=>pdl(1), wavevector=>pdl([0.01]));
+    my $haydock=Photonic::WE::S::Haydock->new(
+	metric=>$ml, nh=>10, keepStates=>1, polarization=>pdl([1])->r2C,
+	epsilon=>$epsilon);
+    my $flo=Photonic::WE::S::Field->new(haydock=>$haydock, nh=>10);
+    my $flv=$flo->rawfield;
+    my $fla=1/$ea;
+    my $flb=1/$eb;
+    #my $fproml=$fla*(1-$gl->f)+$flb*($gl->f);
+    #($fla, $flb)=map {$_/$fproml} ($fla, $flb);
+    my $flx=-4*PI*($fla*(1-$B)+$flb*$B)->transpose;
+    ok(Cagree($flv, $flx), "1D long rawfield");
+}
 
 #View 2D from 1D superlattice. Long wavelength transverse case
-my $Bt=zeroes(1,11)->yvals<5; #2D flat system
-my $epsilont=$ea*(1-$Bt)+$eb*$Bt;
-my $gt=Photonic::Geometry::FromB->new(B=>$Bt); #trans
-my $mt=Photonic::WE::S::Metric->new(geometry=>$gt, epsilon=>pdl(1),
-   wavenumber=>pdl(0.001), wavevector=>pdl([0,0.0001]));
-my $nt=Photonic::WE::S::Haydock->new(
-   metric=>$mt, nh=>10, keepStates=>1, polarization=>pdl([1,0])->r2C,
-   epsilon=>$epsilont);
-my $fto=Photonic::WE::S::Field->new(haydock=>$nt, nh=>10);
-my $ftv=$fto->field;
-my $ftx=r2C(pdl [1, 0]);
-ok(Cagree($ftv, $ftx), "1D trans field");
+{
+    my $Bt=zeroes(1,11)->yvals<5; #2D flat system
+    my $epsilont=$ea*(1-$Bt)+$eb*$Bt;
+    my $gt=Photonic::Geometry::FromB->new(B=>$Bt); #trans
+    my $q=pdl(0.0001);
+    my $k=pdl([0,0.0002]);
+    my $mt=Photonic::WE::S::Metric->new(geometry=>$gt, epsilon=>pdl(1),
+					wavenumber=>$q, wavevector=>$k);
+    my $nt=Photonic::WE::S::Haydock->new(
+	metric=>$mt, nh=>10, keepStates=>1, polarization=>pdl([1,0])->r2C,
+	epsilon=>$epsilont);
+    my $fto=Photonic::WE::S::Field->new(haydock=>$nt, nh=>10);
+    my $ftv=$fto->rawfield;
+    my $f=$gt->f;
+    my $epsM=(1-$f)*$ea+$f*$eb;
+    my $ftx=-4*PI*r2C(pdl [1, 0])->dummy(2,11)*$q**2/($epsM*$q**2-($k**2)->sumover);
+    ok(Cagree($ftv, $ftx), "1D trans rawfield");
+}
