From 5e6d731d530a8af3265b26397fdfe7ed19ce6faa Mon Sep 17 00:00:00 2001 From: Luis Mochan Date: Sun, 22 Sep 2024 11:48:07 -0600 Subject: [PATCH] Add ...ST::Haydock modules and tests --- lib/Photonic/LE/ST/Haydock.pm | 192 +++++++++++++++++++++++++++ lib/Photonic/WE/ST/Haydock.pm | 235 ++++++++++++++++++++++++++++++++++ t/00-load.t | 4 +- t/allh-st.t | 110 ++++++++++++++++ t/allh-west.t | 199 ++++++++++++++++++++++++++++ 5 files changed, 738 insertions(+), 2 deletions(-) create mode 100644 lib/Photonic/LE/ST/Haydock.pm create mode 100644 lib/Photonic/WE/ST/Haydock.pm create mode 100644 t/allh-st.t create mode 100644 t/allh-west.t diff --git a/lib/Photonic/LE/ST/Haydock.pm b/lib/Photonic/LE/ST/Haydock.pm new file mode 100644 index 00000000..6f865959 --- /dev/null +++ b/lib/Photonic/LE/ST/Haydock.pm @@ -0,0 +1,192 @@ +package Photonic::LE::ST::Haydock; +$Photonic::LE::ST::Haydock::VERSION = '0.021'; + +=encoding UTF-8 + +=head1 NAME + +Photonic::LE::S::Haydock + +=head1 VERSION + +version 0.021 + +=head1 COPYRIGHT NOTICE + +Photonic - A perl package for calculations on photonics and +metamaterials. + +Copyright (C) 2016 by W. Luis Mochán + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 1, or (at your option) +any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301 USA + + mochan@fis.unam.mx + + Instituto de Ciencias Físicas, UNAM + Apartado Postal 48-3 + 62251 Cuernavaca, Morelos + México + +=cut + +=head1 SYNOPSIS + + use Photonic::LE::S::Haydock; + my $nr=Photonic::LE::S::Haydock->new(epsilon=>$epsilon, + geometry=>$geometry); + $nr->iterate; + say $nr->iteration; + say $nr->current_a; + say $nr->next_b2; + my $state=$nr->next_state; + +=head1 DESCRIPTION + +Implements calculation of Haydock coefficients and Haydock states for +the calculation of the non retarded dielectric function of arbitrary +periodic N component systems in arbitrary number of dimensions. One +Haydock coefficient at a time. Use k,-k spinors. MQ notes. Assume a +dielectric tensor as input, instead of a dielectric scalar. + +Consumes L, L, +L +- please see those for attributes. + +=head1 ATTRIBUTES SUPPLIED FOR ROLE + +These are provided for roles: + +=over 4 + +=item * geometry GeometryG0 + +A L object defining the geometry of the system, +the characteristic function and the direction of the G=0 vector. Required. + +=item * B ndims dims r G GNorm L scale f + +Accessors handled by geometry (see Photonic::Geometry) + +=item * applyOperator($psi_G) + +Apply the Hamiltonian operator to state in reciprocal space. State is +pm:nx:ny... The operator is the longitudinal dielectric response in +reciprocal-spinor space. + +=item * innerProduct($left, $right) + +Returns the inner (Euclidean) product between states in reciprocal and +spinor space. + +=item * magnitude($psi_G) + +Returns the magnitude of a state by taking the square root of the +inner product of the state with itself. + +=item * changesign + +Returns 0, as there is no need to change sign. + +=item * complexCoeffs + +Haydock coefficients are complex + +=back + +=cut + +use namespace::autoclean; +use PDL::Lite; +use PDL::NiceSlice; +use Carp; +use Photonic::Types -all; +use Photonic::Utils qw(SProd any_complex GtoR RtoG); +use Moo; +use MooX::StrictConstructor; + +has 'geometry'=>(is=>'ro', isa => GeometryG0, + handles=>[qw(B ndims dims r G GNorm L scale f pmGNorm)],required=>1 + ); +has 'complexCoeffs'=>(is=>'ro', init_arg=>undef, default=>1, + documentation=>'Haydock coefficients are complex'); +with 'Photonic::Roles::Haydock', 'Photonic::Roles::UseMask', 'Photonic::Roles::EpsFromGeometry'; + +#Required by Photonic::Roles::Haydock + +sub _build_firstState { #\delta_{G0} + my $self=shift; + my $v=PDL->zeroes(2,@{$self->dims})->r2C; #pm:nx:ny + my $arg=join ',', ':', ("(0)") x $self->ndims; #(0),(0),... ndims+1 times + $v->slice($arg).=1/sqrt(2); + return $v; +} +sub applyOperator { + my $self=shift; + my $psi_G=shift; + my $mask=undef; + $mask=$self->mask if $self->use_mask; + confess "State should be complex" unless any_complex($psi_G); + #Each state is a spinor with two wavefunctions \psi_{k,G} and + #\psi_{-k,G}, thus the index plus or minus k, pm. + #Notation pm=+ or - k, xy=cartesian + #state is pmk:nx:ny... pmGnorm=xy:pmk:nx:ny... + #Multiply by vectors ^G and ^(-G). + #Have to get cartesian out of the way, thread over it and iterate + #over the rest + my $Gpsi_G=($psi_G*$self->pmGNorm->mv(0,-1))->mv(0,-1); #^G |psi> + #the result is complex nx:ny...i:pmk + # Notice that I actually multiply by unit(k-G) instead of + # unit(-k+G) when I use pmGNorm; as I do it two times, the result + # is the same. + #Take inverse Fourier transform over all space dimensions, + #thread over cartesian and pmk indices + #real space ^G|psi> + my $Gpsi_R=GtoR($Gpsi_G, $self->ndims, 0); # $Gpsi_R is nx:ny...i:pmk + # $self->epsilon is j:i:nx:ny... + #Multiply by the dielectric function in Real Space. Thread + #cartesian and pm indices + my $eGpsi_R= ( $self->epsilon # Epsilon is tensorial j:i:nx:ny... + * $Gpsi_R->mv(-2,0)->dummy(1) # j:i:nx:ny...:pmk + )->sumover #matrix product i:nx:ny...:pmk + ->mv(0,-2); + #$eGpsi_R is nx:ny,...i:pmk + #Transform to reciprocal space + my $eGpsi_G=RtoG($eGpsi_R, $self->ndims, 0)->mv(-1,0); + #$eGpsi_G is pmk:nx:ny...:i + #Scalar product with pmGnorm: i:pm:nx:ny... + my $GeGpsi_G=($eGpsi_G*$self->pmGNorm->mv(0,-1)) #^Ge^G|psi> + # pmk:nx:ny...:i + # Move cartesian to front and sum over + ->mv(-1,0)->sumover; #^G.epsilon^G|psi> + # $GeGpsi_G is pm:nx:ny. $mask=nx:ny + $GeGpsi_G=$GeGpsi_G*$mask->(*1) if defined $mask; + return $GeGpsi_G; +} +sub innerProduct { + #ignore self + return SProd($_[1], $_[2]); +} +sub magnitude { + my $self=shift; + my $psi=shift; + return $self->innerProduct($psi, $psi)->abs->sqrt; +} +sub changesign { #don't change sign + return 0; +} + +__PACKAGE__->meta->make_immutable; + +1; diff --git a/lib/Photonic/WE/ST/Haydock.pm b/lib/Photonic/WE/ST/Haydock.pm new file mode 100644 index 00000000..d3f1bc2d --- /dev/null +++ b/lib/Photonic/WE/ST/Haydock.pm @@ -0,0 +1,235 @@ +package Photonic::WE::ST::Haydock; +$Photonic::WE::ST::Haydock::VERSION = '0.021'; + +=encoding UTF-8 + +=head1 NAME + +Photonic::WE::ST::Haydock + +=head1 VERSION + +version 0.021 + +=head1 SYNOPSIS + + use Photonic::WE::ST::Haydock; + my $nr=Photonic::WE::ST::Haydock->new(metric=>$g, polarization=>$p, nh=>$nh); + $nr->iterate; + say $nr->iteration; + say $nr->current_a; + say $nr->next_b2; + my $state=$nr->next_state; + +=head1 DESCRIPTION + +Implements calculation of Haydock coefficients and Haydock states for +the calculation of the retarded dielectric function of arbitrary +periodic systems in arbitrary number of dimensions, one +Haydock coefficient at a time. It uses the wave equation and the spinor representation. + +Consumes L, L, +L +- please see those for attributes. + +=head1 ATTRIBUTES + +=over 4 + +=item * metric Photonic::WE::ST::Metric + +A L object defining the geometry of the +system, the characteristic function, the wavenumber, wavevector and +host dielectric function. Required in the initializer. + +=item * B ndims dims epsilon + +Accessors handled by metric (see L) + +=item * polarization complex PDL + +A non null vector defining the complex direction of the macroscopic +field. + +=item * normalizedPolarization + +The polarisation, normalised + +=back + +=head1 METHODS + +=over 4 + +=item * applyMetric($psi) + +Returns the result of applying the metric to the state $psi. + +=back + +=head1 ATTRIBUTES SUPPLIED FOR ROLE + +=over 4 + +=item * applyOperator($psi_G) + +Apply the 'Hamiltonian' operator to state $psi_G. State is +ri,xy,pn,nx,ny... The Hamiltonian is the metric followed by the +dielectric esponse relative to the reference response. + +=item * innerProduct($left, $right) + +Returns the inner Euclidean product between states with the metric. + +=item * magnitude($psi) + +Returns the magnitude of a state as the square root of +the magnitude of inner product of the state with itself. + +=item * changesign + +Returns 0, as there is no need to change sign. + +=item * complexCoeffs + +Haydock coefficients are complex + +=back + +=cut + +=head1 COPYRIGHT NOTICE + +Photonic - A perl package for calculations on photonics and +metamaterials. + +Copyright (C) 2016 by W. Luis Mochán + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 1, or (at your option) +any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301 USA + + mochan@fis.unam.mx + + Instituto de Ciencias Físicas, UNAM + Apartado Postal 48-3 + 62251 Cuernavaca, Morelos + México + +=cut + +use v5.36; +use namespace::autoclean; +use PDL::Lite; +use PDL::NiceSlice; +use Carp; +use Photonic::Types -all; +use Photonic::Utils qw(VSProd any_complex GtoR RtoG); +use Moo; +use MooX::StrictConstructor; + +has 'metric'=>(is=>'ro', isa => InstanceOf['Photonic::WE::ST::Metric'], + handles=>{B=>'B', ndims=>'ndims', dims=>'dims', + geometry=>'geometry', epsilonR=>'epsilon'}, + required=>1); +has 'polarization' =>(is=>'ro', required=>1, isa=>PDLComplex); +has 'normalizedPolarization' =>(is=>'ro', isa=>PDLComplex, + init_arg=>undef, writer=>'_normalizedPolarization'); +has 'complexCoeffs'=>(is=>'ro', init_arg=>undef, default=>1, + documentation=>'Haydock coefficients are complex'); +with 'Photonic::Roles::Haydock', 'Photonic::Roles::UseMask', 'Photonic::Roles::EpsFromGeometry'; + +sub applyOperator { + my $self=shift; + my $psi=shift; #psi is xy:pm:nx:ny... + my $mask=undef; + $mask=$self->mask if $self->use_mask; + my $id=PDL::MatrixOps::identity($self->ndims); + my $gpsi=$self->applyMetric($psi); + # gpsi is xy:pm:nx:ny. Get cartesian and pm out of the way and + my $gpsi_r=GtoR($gpsi, $self->ndims, 2)->mv(1,-1); #xy:nx:ny:pm + #my $gpsi_r=GtoR($gpsi, $self->ndims, 2)->mv(0,-1)->mv(0,-1); + ##nx:ny:xy:pm + my $H=($self->epsilonR*$id-$self->epsilon)/$self->epsilonR; #j:i:nx:ny (j:i cartesian) + my $Hgpsi_r=( + $H #j:i:nx:ny + *$gpsi_r->(,*1) #j:i:nx:ny:pm + )->sumover; #i:nx:ny:pm +# ->mv(0,-2); #nx:ny:i:pm QUITAR + #Transform to reciprocal space, move xy and pm back and make complex, + my $psi_G=RtoG( + $Hgpsi_r #i:nx:ny:pm + ->mv(-1,1), #i:pm:nx:ny + # ->mv(-1,0), # QUITAR + $self->ndims, 2); + #Apply mask + #psi_G is xy:pm:nx:ny mask is nx:ny + $psi_G=$psi_G*$mask->(*1,*1) if defined $mask; #use dummies for xy:pm + return $psi_G; +} + +sub applyMetric { + my $self=shift; + my $psi=shift; + #psi is xy:pm:nx:ny + my $g=$self->metric->value; + #$g is xy:xy:pm:nx:ny + #real or complex matrix times complex vector + my $gpsi=($g # j:i:pm:nx:nx + *$psi(:,*1)) #j:i:pm:nx:ny + ->sumover; #i:pm:nx:ny + return $gpsi; +} + +sub innerProduct { #Return Symmetric product with metric + my $self=shift; + my $psi1=shift; + my $psi2=shift; + my $gpsi2=$self->applyMetric($psi2); + return VSProd($psi1, $gpsi2); +} + + +sub magnitude { + my $self=shift; + my $psi=shift; + return $self->innerProduct($psi, $psi)->abs->sqrt; +} +sub changesign { #don't change sign + return 0; +} + +sub _build_firstState { #\delta_{G0} + my $self=shift; + my $v=PDL->zeroes(@{$self->dims})->r2C; #nx:ny... + my $arg=join ',', ("(0)") x $self->ndims; #(0),(0),... ndims times + $v->slice($arg).=1/sqrt(2); + my $e=$self->polarization; #xy + my $d=$e->dim(0); + confess "Polarization has wrong dimensions. " . + " Should be $d-dimensional complex vector, got ($e)." + unless any_complex($e) && $e->dim(0)==$d; + my $modulus2=$e->abs2->sumover; + confess "Polarization should be non null" unless + $modulus2 > 0; + $e=$e/sqrt($modulus2); + $self->_normalizedPolarization($e); + #Use same helicity for k and for -k, conjugate polarization, for allowing chirality + my $phi=pdl($e, $e->conj)*$v->(*1,*1); #initial state ordinarily normalized + # xy:pm:nx:ny + return $phi; +} + +__PACKAGE__->meta->make_immutable; + +1; diff --git a/t/00-load.t b/t/00-load.t index bc39e93e..924da49a 100644 --- a/t/00-load.t +++ b/t/00-load.t @@ -56,6 +56,7 @@ Photonic::LE::S::EpsL Photonic::LE::S::EpsTensor Photonic::LE::S::Field Photonic::LE::S::Haydock +Photonic::LE::ST::Haydock Photonic::Roles::EpsFromGeometry Photonic::Roles::EpsL Photonic::Roles::EpsTensor @@ -78,15 +79,14 @@ Photonic::WE::S::Green Photonic::WE::S::GreenP Photonic::WE::S::Haydock Photonic::WE::S::Metric +Photonic::WE::ST::Haydock ); #Photonic::LE::ST::EpsL #Photonic::LE::ST::EpsTensor #Photonic::LE::ST::Field -#Photonic::LE::ST::Haydock #Photonic::WE::ST::Field #Photonic::WE::ST::Green #Photonic::WE::ST::GreenP -#Photonic::WE::ST::Haydock #Photonic::WE::ST::Metric diff --git a/t/allh-st.t b/t/allh-st.t new file mode 100644 index 00000000..c1c6122e --- /dev/null +++ b/t/allh-st.t @@ -0,0 +1,110 @@ +=head1 COPYRIGHT NOTICE + +Photonic - A perl package for calculations on photonics and +metamaterials. + +Copyright (C) 2016 by W. Luis Mochán + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 1, or (at your option) +any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301 USA + + mochan@fis.unam.mx + + Instituto de Ciencias Físicas, UNAM + Apartado Postal 48-3 + 62251 Cuernavaca, Morelos + México + +=cut + +use strict; +use warnings; +use PDL; +use Photonic::Geometry::FromEpsilonTensor; +use Photonic::LE::ST::Haydock; + +use Test::More; +use lib 't/lib'; +use TestUtils; + +my $fn = make_fn(); + +#Check haydock coefficients for simple 1D system +my ($ea, $eb)=(1+2*i, 3+4*i); +my $f=6/11; +my $eps=r2C($ea*(zeroes(1,1,11)->zvals<5)+ $eb*(zeroes(1,1,11)->zvals>=5)); +my $g=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$eps, Direction0=>pdl([1])); +my $a=Photonic::LE::ST::Haydock->new(geometry=>$g, nh=>10); +$a->run; +my $as=$a->as; +my $bs=$a->bs; +my $b2s=$a->b2s; +is($a->iteration, 2, "Number of iterations 1D longitudinal"); +ok(Cagree($b2s->slice("(0)"), r2C(1)), "1D L b_0^2"); +ok(Cagree($as, pdl([$ea*(1-$f)+$eb*$f, $ea*$f+$eb*(1-$f)])), "1D L a"); +ok(Cagree($b2s->slice("(1)"), ($eb-$ea)**2*$f*(1-$f)), "1D L b_1^2"); +ok(Cagree($b2s, $bs**2), "1D L b2==b^2"); + +#View 1D system as 2D, isotropic but tensorial. Transverse direction +my $epst=($ea*identity(2)*(zeroes(2,2,11,1)->zvals<5) + +$eb*identity(2)*(zeroes(2,2,11,1)->zvals>=5) + )->mv(2,3); +my $gt=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$epst, Direction0=>pdl([1,0])); #trans +my $at=Photonic::LE::ST::Haydock->new(geometry=>$gt, nh=>10); +$at->run; +my $ast=$at->as; +my $bst=$at->bs; +my $b2st=$at->b2s; +is($at->iteration, 1, "Number of iterations 1D trans"); +ok(Cagree($b2st->slice("(0)"), 1), "1D T b_0^2"); +ok(Cagree($ast->slice("(0)"), $ea*(1-$f)+$eb*$f), "1D T a_0"); +ok(Cagree($b2st, $bst**2), "1D T b2==b^2"); + +{ + #check reorthogonalize with square array + my $epss=($eb*(zeroes(15,15)->rvals<5)+$ea*(zeroes(15,15)->rvals>=5))->slice("*2,*2") + *identity(2); + my $gs=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$epss, Direction0=>pdl([1,0]), L=>pdl(1,1)); + my $als=Photonic::LE::ST::Haydock + ->new(geometry=>$gs, nh=>2*15*15, reorthogonalize=>1, + accuracy=>machine_epsilon(), noise=>3*machine_epsilon(), + normOp=>$eb->abs->sumover->sumover # sum of abs of matrix elements + ); + $als->run; + ok($als->iteration <= 15*15, + "No more iterations than dimensions. Square. States in mem."); + diag("Actual iterations: " .$als->iteration + . " Actual orthogonalizations: " . $als->orthogonalizations); +} +{ + #check reorthogonalize with square array. Data in file. + my $epss=($eb*(zeroes(15,15)->rvals<5)+$ea*(zeroes(15,15)->rvals>=5))->slice("*2,*2") + *identity(2); + my $gs=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$epss, Direction0=>pdl([1,0]), L=>pdl(1,1)); + my $als=Photonic::LE::ST::Haydock + ->new(geometry=>$gs, nh=>2*15*15, reorthogonalize=>1, + accuracy=>machine_epsilon(), noise=>3*machine_epsilon(), + normOp=>$eb->abs->sumover->sumover, stateFN=>$fn); + $als->run; + ok($als->iteration <= 15*15, + "No more iterations than dimensions. Square. States in file"); + diag("Actual iterations: " .$als->iteration + . " Actual orthogonalizations: " . $als->orthogonalizations); +} + +done_testing; diff --git a/t/allh-west.t b/t/allh-west.t new file mode 100644 index 00000000..9e094dbf --- /dev/null +++ b/t/allh-west.t @@ -0,0 +1,199 @@ +=head1 COPYRIGHT NOTICE + +Photonic - A perl package for calculations on photonics and +metamaterials. + +Copyright (C) 2016 by W. Luis Mochán + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 1, or (at your option) +any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301 USA + + mochan@fis.unam.mx + + Instituto de Ciencias Físicas, UNAM + Apartado Postal 48-3 + 62251 Cuernavaca, Morelos + México + +=cut + +use strict; +use warnings; +use PDL; +use Photonic::Geometry::FromEpsilonTensor; +use Photonic::WE::ST::Metric; +use Photonic::WE::ST::Haydock; + +use Test::More; +use lib 't/lib'; +use TestUtils; + +my $fn = make_fn(); + +#Check haydock coefficients for simple 1D system +my ($ea, $eb)=(1+2*i, 3+4*i); +my $f=6/11; +my $eps=$ea*(zeroes(1,1,11)->zvals<5)+ $eb*(zeroes(1,1,11)->zvals>=5); +my $g=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$eps); +my $m=Photonic::WE::ST::Metric->new( + geometry=>$g, epsilon=>pdl(1), wavenumber=>pdl(2), wavevector=>pdl([1]) + ); +my $a=Photonic::WE::ST::Haydock->new(metric=>$m, + polarization=>pdl([1])->r2C, nh=>10); +$a->run; +my $as=$a->as; +my $bs=$a->bs; +my $b2s=$a->b2s; +is($a->iteration, 2, "Number of iterations 1D longitudinal x"); +ok(Cagree($b2s, pdl([1, ($eb-$ea)**2*$f*(1-$f)])), "1D L b^2"); +ok(Cagree($as, pdl([(1-$ea)*(1-$f)+(1-$eb)*$f, (1-$ea)*$f+(1-$eb)*(1-$f)])), "1D L a"); +ok(Cagree($b2s, $bs**2), "1D L b2==b^2"); + +#Check haydock coefficients for simple 1D system other longitudinal y +my $eps1l=($ea*(zeroes(11,1)->xvals<5)+ $eb*(zeroes(11,1)->xvals>=5))->slice("*1,*1") + *identity(2); +my $g1l=Photonic::Geometry::FromEpsilonTensor + ->new(epsilon=>$eps1l); +my $m1l=Photonic::WE::ST::Metric->new( + geometry=>$g1l, epsilon=>pdl(1), wavenumber=>pdl(.0002), + wavevector=>pdl([0,.000001]) + ); +my $a1l=Photonic::WE::ST::Haydock->new(metric=>$m1l, + polarization=>pdl([0,1])->r2C, nh=>10, smallH=>1e-4); +$a1l->run; +my $as1l=$a1l->as; +my $b2s1l=$a1l->b2s; +is($a1l->iteration, 1, "Number of iterations 1D long y"); +ok(Cagree(($b2s1l->slice("(0)")), 1), "1D L b_0^2"); +ok(Cagree($as1l->slice("(0)"), (1-$ea)*(1-$f)+(1-$eb)*$f), "1D L a_0"); + +#Check haydock coefficients for simple 1D system transverse prop x pol y +my $epst=($ea*(zeroes(11,1)->xvals<5)+ $eb*(zeroes(11,1)->xvals>=5))->slice("*1,*1") + *identity(2); +my $gt=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$epst); +my $mt=Photonic::WE::ST::Metric->new( + geometry=>$gt, epsilon=>pdl(1), wavenumber=>pdl(.0002), + wavevector=>pdl([.000001,0]) + ); +my $at=Photonic::WE::ST::Haydock->new(metric=>$mt, + polarization=>pdl([0,1])->r2C, nh=>10, smallH=>1e-4); +$at->run; +my $ast=$at->as; +my $bst=$at->bs; +my $b2st=$at->b2s; +is($at->iteration, 1, "Number of iterations 1D trans"); +ok(Cagree($b2st->slice("(0)"), 1), "1D L b_0^2"); +ok(Cagree($ast->slice("(0)"), (1-$ea)*(1-$f)+(1-$eb)*$f), "1D L a_0"); + +#check reorthogonalize with square array +my $Bs=zeroes(9,9)->rvals<4; +my $epss=r2C((1-$Bs)+5*$Bs)->slice("*1,*1")*identity(2); +my $gs=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$epss, L=>pdl(1,1)); +my $ms=Photonic::WE::ST::Metric->new(geometry=>$gs, epsilon=>pdl(1), + wavenumber=>pdl(.01), wavevector=>pdl([.001,0])); +my $als=Photonic::WE::ST::Haydock + ->new(metric=>$ms, polarization=>r2C(pdl([0,1])), nh=>2*9*9, + reorthogonalize=>1, accuracy=>machine_epsilon(), + noise=>1e1*machine_epsilon(), normOp=>1e0, smallH=>1e-7); +$als->run; +ok($als->iteration < 2*9*9, + "No more iterations than dimensions. Square. Long wavelength."); +diag("Actual iterations: " . $als->iteration + . " Actual orthogonalizations: ", $als->orthogonalizations); + +{ + #check reorthogonalize with even numbers + my $Be=zeroes(10,10)->rvals<4; + my $epse=r2C((1-$Be)+5*$Be)->slice("*1,*1")*identity(2); + my $ge=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$epse, L=>pdl(1,1)); + my $me=Photonic::WE::ST::Metric->new( + geometry=>$ge, epsilon=>pdl(1), wavenumber=>pdl(.01), + wavevector=>pdl([.001,0])); + my $ale=Photonic::WE::ST::Haydock + ->new(metric=>$me, polarization=>r2C(pdl([0,1])), nh=>2*15*15, + reorthogonalize=>1, accuracy=>machine_epsilon(), + noise=>1e1*machine_epsilon(), normOp=>1e0, smallH=>1e-7, + use_mask=>1); + $ale->run; + ok($ale->iteration < 2*10*10, + "No more iterations than dimensions. Square. Long wavelength. Even."); + diag("Actual iterations: " . $ale->iteration + . " Actual orthogonalizations: ", $ale->orthogonalizations); +} + +{ + #check reorthogonalize with even numbers + my $Be=zeroes(10,10)->rvals<4; + my $epse=r2C((1-$Be)+5*$Be)->slice("*1,*1")*identity(2); + my $ge=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$epse, L=>pdl(1,1)); + my $me=Photonic::WE::ST::Metric->new( + geometry=>$ge, epsilon=>pdl(1), wavenumber=>pdl(3.6), + wavevector=>pdl([1.01*3.6,0])); + my $ale=Photonic::WE::ST::Haydock + ->new(metric=>$me, polarization=>r2C(pdl([0,1])), nh=>2*15*15, + reorthogonalize=>1, accuracy=>machine_epsilon(), + noise=>1e1*machine_epsilon(), normOp=>1e0, smallH=>1e-7, + use_mask=>1); + $ale->run; + ok($ale->iteration < 2*10*10, + "No more iterations than dimensions. Square. Short wavelength. Even."); + diag("Actual iterations: " . $ale->iteration + . " Actual orthogonalizations: ", $ale->orthogonalizations); +} + +{ + #check reorthogonalize with even numbers + my $Be=zeroes(10,10)->rvals<4; + my $epse=r2C((1-$Be)+5*$Be)->slice("*1,*1")*identity(2); + my $ge=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$epse, L=>pdl(1,1)); + my $me=Photonic::WE::ST::Metric->new( + geometry=>$ge, epsilon=>pdl(1), wavenumber=>pdl(3.6), + wavevector=>pdl([1.01*3.6,0])); + my $ale=Photonic::WE::ST::Haydock + ->new(metric=>$me, polarization=>r2C(pdl([0,1])), nh=>2*15*15, + reorthogonalize=>1, accuracy=>machine_epsilon(), + noise=>1e1*machine_epsilon(), normOp=>1e0, smallH=>1e-7, + use_mask=>1, stateFN=>$fn); + $ale->run; + ok($ale->iteration < 2*10*10, + "No more iterations than dimensions. Square. Short wavelength. Even. Disk"); + diag("Actual iterations: " . $ale->iteration + . " Actual orthogonalizations: ", $ale->orthogonalizations); +} + +done_testing; + +__END__ + +#check reorthogonalize again with square array +my $B1s=zeroes(21,21)->rvals>2; #21,21, 2 +my $eps1s=r2C(10*(1-$B1s)+1*$B1s); +my $g1s=Photonic::Geometry::FromEpsilon->new(epsilon=>$eps1s, L=>pdl(1,1)); +my $m1s=Photonic::WE::S::Metric->new(geometry=>$g1s, epsilon=>pdl(10), + wavenumber=>pdl(3.6), wavevector=>pdl([1.01*3.6,0])); +my $al1s=Photonic::WE::S::Haydock + ->new(metric=>$m1s, polarization=>r2C(pdl([0,1])), nh=>3*21*21, + reorthogonalize=>1, accuracy=>machine_epsilon(), + noise=>1e3*machine_epsilon(), normOp=>1e0, smallH=>1e-7); +$al1s->run; +ok($al1s->iteration <= 2*21*21, "No more iterations than dimensions"); +diag("Actual iterations: " . $al1s->iteration + . " Actual orthogonalizations: ", $al1s->orthogonalizations); + + +#foreach(@$st){ +# my $pr=$als->innerProduct($_, $st->[0]); +# print "$pr\n"; +#}