diff --git a/MANIFEST b/MANIFEST
index cb171071..ea4bf2d1 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -57,9 +57,9 @@ lib/Photonic/WE/S/Green.pm
lib/Photonic/WE/S/GreenP.pm
lib/Photonic/WE/S/Haydock.pm
lib/Photonic/WE/S/Metric.pm
-# lib/Photonic/WE/ST/Field.pm
-# lib/Photonic/WE/ST/Green.pm
-# lib/Photonic/WE/ST/GreenP.pm
+lib/Photonic/WE/ST/Field.pm
+lib/Photonic/WE/ST/Green.pm
+lib/Photonic/WE/ST/GreenP.pm
lib/Photonic/WE/ST/Haydock.pm
lib/Photonic/WE/ST/Metric.pm
LICENSE
@@ -89,13 +89,13 @@ t/epst-nr2.t
t/epst-st.t
t/epst-s.t
t/epst-wer2.t
-# t/epst-west.t
+t/epst-west.t
t/epst-wes.t
t/field-lenr2.t
t/field-lest.t
t/field-les.t
t/field-wer2.t
-# t/field-west.t
+t/field-west.t
t/field-wes.t
t/geometry.t
t/greenp.t
diff --git a/lib/Photonic/WE/ST/Field.pm b/lib/Photonic/WE/ST/Field.pm
new file mode 100644
index 00000000..52119965
--- /dev/null
+++ b/lib/Photonic/WE/ST/Field.pm
@@ -0,0 +1,117 @@
+package Photonic::WE::ST::Field;
+$Photonic::WE::ST::Field::VERSION = '0.022';
+
+=encoding UTF-8
+
+=head1 NAME
+
+Photonic::WE::ST::Field
+
+=head1 VERSION
+
+version 0.022
+
+=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::WE::ST::Field;
+ my $nrf=Photonic::WE::ST::Field->new(...);
+ my $field=$nrf->field;
+
+=head1 DESCRIPTION
+
+Calculates the non retarded electric field for a given fixed
+Photonic::Geometry structure and given dielectric functions of
+the components.
+
+Consumes L<Photonic::Roles::Field>
+- please see for attributes.
+
+=cut
+
+use namespace::autoclean;
+use PDL::Lite;
+use PDL::NiceSlice;
+use Photonic::WE::ST::Haydock;
+use Photonic::Utils qw(cgtsv GtoR linearCombineIt);
+use Photonic::Types -all;
+use Moo;
+use MooX::StrictConstructor;
+
+with 'Photonic::Roles::Field';
+
+sub BUILD {
+ my $self=shift;
+ $self->haydock->run unless $self->haydock->iteration;
+}
+
+sub _build_field {
+ my $self=shift;
+ my $as=$self->haydock->as;
+ my $bs=$self->haydock->bs;
+ my $cs=$self->haydock->cs;
+ my $nh=$self->nh; #desired number of Haydock terms
+ #don't go beyond available values.
+ $nh=$self->haydock->iteration if $nh>$self->haydock->iteration;
+ # calculate using lapack for tridiag system
+ my $diag = 1-$as->(0:$nh-1);
+ # rotate complex zero from first to last element.
+ my $subdiag = -$bs->(0:$nh-1)->rotate(-1);
+ my $supradiag =-$cs->(0:$nh-1)->rotate(-1);
+ my $rhs=PDL->zeroes($nh); #build a nh pdl
+ $rhs->slice((0)).=1;
+ $rhs=$rhs->r2C;
+ #coefficients of g^{-1}E
+ my $giEs= cgtsv($subdiag, $diag, $supradiag, $rhs);
+ #states are xy,pm,nx,ny...
+ my $stateit=$self->haydock->states;
+ my $ndims=$self->haydock->B->ndims; # num. of dims of space
+ #field is xy,pm,nx,ny...
+ my $field_G=linearCombineIt($giEs, $stateit); #En ^G|psi_n>
+ my $Es=$self->haydock->applyMetric($field_G);
+ #Comment as normalization below makes it useless
+ #$Es*=$bs->((0))/$self->haydock->metric->epsilon;
+ my $Esp=$Es(:,(0)); # choose +k spinor component.
+ my $e_0=1/($Esp->slice(":" . ",(0)" x $ndims)
+ *$self->haydock->polarization->conj)->sumover;
+ # Normalize result so macroscopic field is 1.
+ $Esp*=$e_0;
+ $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 to have unit macroscopic field
+ return $field_R; #result is xy,nx,ny...
+}
+
+__PACKAGE__->meta->make_immutable;
+
+1;
diff --git a/lib/Photonic/WE/ST/Green.pm b/lib/Photonic/WE/ST/Green.pm
new file mode 100644
index 00000000..44940312
--- /dev/null
+++ b/lib/Photonic/WE/ST/Green.pm
@@ -0,0 +1,246 @@
+package Photonic::WE::ST::Green;
+$Photonic::WE::ST::Green::VERSION = '0.022';
+
+=encoding UTF-8
+
+=head1 NAME
+
+Photonic::WE::ST::Green
+
+=head1 VERSION
+
+version 0.022
+
+=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::WE::ST::Green;
+ my $G=Photonic::WE::ST::Green->new(metric=>$m, nh=>$nh);
+ my $GreenTensor=$G->greenTensor;
+ my $WaveTensor=$G->waveOperator;
+ my $EpsTensor=$G->epsilonTensor;
+
+=head1 DESCRIPTION
+
+Calculates the retarded green's tensor for a given fixed
+Photonic::WE::ST::Metric structure as a function of the dielectric
+functions of the components.
+
+=head1 ATTRIBUTES
+
+=over 4
+
+=item * keepStates
+
+Value of flag to keep Haydock states in Haydock calculations (default 0)
+
+=item * metric
+
+L<Photonic::WE::ST::Metric> describing the structure and some parameters.
+
+=item * nh
+
+The maximum number of Haydock coefficients to use.
+
+=item * smallH, smallE
+
+Criteria of convergence of Haydock coefficients and continued
+fraction. 0 means don't check. (default 1e-7)
+
+=item * haydock
+
+Array of L<Photonic::WE::ST::Haydock> structures, one for each polarization
+
+=item * reorthogonalize
+
+Reorthogonalize haydock flag
+
+=item * greenP
+
+Array of L<Photonic::WE::ST::GreenP> structures, one for each direction.
+
+=item * greenTensor
+
+The Green's tensor calculated
+
+=item * nhActual
+
+The actual number of Haydock coefficients used in the last calculation
+
+=item * converged
+
+Flags that the last calculation converged before using up all coefficients
+
+=item * waveOperator
+
+The macroscopic wave operator of the last operation
+
+=item * epsilonTensor
+
+The macroscopic dielectric tensor
+
+=back
+
+=cut
+
+use namespace::autoclean;
+use PDL::Lite;
+use PDL::NiceSlice;
+use Photonic::WE::ST::Haydock;
+use Photonic::WE::ST::GreenP;
+use Photonic::Types -all;
+use Photonic::Utils qw(tensor make_haydock make_greenp wave_operator any_complex triangle_coords);
+use List::Util qw(all);
+use Moo;
+use MooX::StrictConstructor;
+
+has 'nh' =>(is=>'ro', isa=>Num, required=>1,
+ documentation=>'Desired no. of Haydock coefficients');
+has 'smallH'=>(is=>'ro', isa=>Num, required=>1, default=>1e-7,
+ documentation=>'Convergence criterium for Haydock coefficients');
+has 'smallE'=>(is=>'ro', isa=>Num, required=>1, default=>1e-7,
+ documentation=>'Convergence criterium for use of Haydock coeff.');
+has 'metric'=>(is=>'ro', isa => InstanceOf['Photonic::WE::ST::Metric'],
+ handles=>[qw(geometry ndims dims)],required=>1);
+
+has 'haydock' =>(is=>'lazy', isa=>ArrayRef[Haydock],
+ init_arg=>undef,
+ documentation=>'Array of Haydock calculators');
+has 'greenP'=>(is=>'lazy', isa=>ArrayRef[InstanceOf['Photonic::WE::ST::GreenP']],
+ init_arg=>undef,
+ documentation=>'Array of projected G calculators');
+has 'converged'=>(is=>'ro', init_arg=>undef, writer=>'_converged',
+ documentation=>
+ 'All greenP evaluations converged');
+has 'greenTensor'=>(is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'Greens Tensor');
+has 'reorthogonalize'=>(is=>'ro', required=>1, default=>0,
+ documentation=>'Reorthogonalize haydock flag');
+has 'waveOperator' => (is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'Wave operator');
+has 'epsilonTensor' => (is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'macroscopic response');
+
+with 'Photonic::Roles::KeepStates', 'Photonic::Roles::UseMask';
+
+has 'cHaydock' =>(
+ is=>'lazy', isa=>ArrayRef[Haydock],
+ init_arg=>undef,
+ documentation=>'Array of Haydock calculators for complex projection');
+
+has 'cGreenP'=>(
+ is=>'lazy', isa=>ArrayRef[InstanceOf['Photonic::WE::ST::GreenP']],
+ init_arg=>undef,
+ documentation=>'Array of projected G calculators for complex projection');
+
+has 'symmetric' => (
+ is=>'ro', required=>1, default=>0,
+ documentation=>'Flags only symmetric part required');
+
+sub _build_greenTensor {
+ my $self=shift;
+ my $symmetric=tensor(
+ pdl([map $_->Gpp, @{$self->greenP}]),
+ $self->geometry->unitDyadsLU,
+ my $nd=$self->geometry->ndims, 2);
+ $self->_converged(all { $_->converged } @{$self->greenP});
+ return $symmetric if $self->symmetric;
+ my $greenPc = pdl map $_->Gpp,
+ my @cGP=@{$self->cGreenP}; #Green's projections along complex directions.
+ $self->_converged(all { $_->converged } $self, @cGP);
+ my $asy=$symmetric->zeroes; #xy,xy, $ndx$nd
+ my $cpairs=$self->geometry->cUnitPairs->mv(1,-1);
+ my $indexes = triangle_coords($nd);
+ $indexes = $indexes->mv(-1,0)->whereND( ($indexes(0) <= $nd-2)->((0)) )->mv(0,-1); # first index only up to $nd-2, mv because whereND takes dims off bottom
+ $asy->indexND($indexes) .= #$asy is xy,xy. First index is column
+ $greenPc-
+ ($cpairs->conj->(*1) # column, row
+ *$cpairs->(:,*1)
+ *$symmetric)->sumover->sumover
+ ;
+ $asy *= PDL->i();
+ # This is wrong: $asy -= $asy->transpose;
+ $asy = $asy-$asy->transpose;
+ $symmetric+$asy;
+}
+
+
+sub _build_haydock { # One Haydock coefficients calculator per direction0
+ my ($self) = @_;
+ make_haydock($self, 'Photonic::WE::ST::Haydock', $self->geometry->unitPairs, 0, qw(reorthogonalize use_mask mask));
+}
+
+sub _build_cHaydock {
+ # One Haydock coefficients calculator per complex polarization
+ my $self=shift;
+ make_haydock($self, 'Photonic::WE::ST::Haydock', $self->geometry->cUnitPairs, 0,
+ qw(reorthogonalize use_mask mask));
+}
+
+sub _build_greenP {
+ make_greenp(shift, 'Photonic::WE::ST::GreenP');
+}
+
+sub _build_cGreenP {
+ make_greenp(shift, 'Photonic::WE::ST::GreenP', 'cHaydock');
+}
+
+sub _build_waveOperator {
+ my $self=shift;
+ wave_operator($self->greenTensor, $self->geometry->ndims);
+}
+
+sub _build_epsilonTensor {
+ my $self=shift;
+ my $wave=$self->waveOperator;
+ my $q=$self->metric->wavenumber;
+ my $q2=$q*$q;
+ my $k=$self->metric->wavevector;
+ my ($k2, $kk);
+ if(any_complex($q, $k)){
+ #Make both complex
+ $_ = PDL::r2C($_) for $q, $k;
+ $k2=($k*$k)->sumover; #inner
+ $kk=$k->(:,*1)*$k->(*1); #outer
+ } else {
+ $k2=$k->inner($k);
+ $kk=$k->outer($k);
+ }
+ my $id=PDL::MatrixOps::identity($k);
+ $wave+$k2/$q2*$id - $kk/$q2;
+}
+
+__PACKAGE__->meta->make_immutable;
+
+1;
+
+__END__
diff --git a/lib/Photonic/WE/ST/GreenP.pm b/lib/Photonic/WE/ST/GreenP.pm
new file mode 100644
index 00000000..1bcb1dc4
--- /dev/null
+++ b/lib/Photonic/WE/ST/GreenP.pm
@@ -0,0 +1,173 @@
+package Photonic::WE::ST::GreenP;
+$Photonic::WE::ST::GreenP::VERSION = '0.022';
+
+=encoding UTF-8
+
+=head1 NAME
+
+Photonic::WE::ST::GreenP
+
+=head1 VERSION
+
+version 0.022
+
+=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::WE::ST::GreenP;
+ my $green=Photonic::WE::ST::GreenP->new(haydock=>$h, nh=>$nh);
+ my $greenProjection=$green->Gpp;
+ my $WaveProjection=$green->waveOperator;
+ my $EpsTensor=$green->epsilon;
+
+=head1 DESCRIPTION
+
+Calculates the dielectric function for a given fixed
+L<Photonic::WE::ST::Haydock> structure as a function of the dielectric
+functions of the components.
+
+=head1 ATTRIBUTES
+
+=over 4
+
+=item * haydock
+
+The L<Photonic::WE::ST::Haydock> structure (required).
+
+=item * nh
+
+The maximum number of Haydock coefficients to use.
+
+=item * smallE
+
+Criteria of convergence. 0 means don't check. (defaults to 1e-7)
+
+=item * u
+
+The spectral variable used in the calculation
+
+=item * nhActual
+
+The actual number of Haydock coefficients used in the calculation
+
+=item * converged
+
+Flags that the calculation converged before using up all coefficients
+
+=item * waveOperator
+
+The macroscopic wave operator calculated from the metric.
+
+NOTE: Only works along principal directions, as it treats Green's
+function as scalar.
+
+=item * epsilon
+
+The macroscopic dielectric projection
+
+NOTE: Only works for polarizations along principal directions.
+
+=back
+
+=cut
+
+use namespace::autoclean;
+use PDL::Lite;
+use Photonic::WE::ST::Haydock;
+use Photonic::Types -all;
+use Photonic::Utils qw(lentzCF);
+use List::Util qw(min);
+use Moo;
+use MooX::StrictConstructor;
+
+has 'nh' =>(is=>'ro', isa=>Num, required=>1,
+ documentation=>'Desired no. of Haydock coefficients');
+has 'smallE'=>(is=>'ro', isa=>Num, required=>1, default=>1e-7,
+ documentation=>'Convergence criterium for use of Haydock coeff.');
+has 'haydock' =>(is=>'ro', isa=>Haydock, required=>1);
+has 'nhActual'=>(is=>'ro', isa=>Num, init_arg=>undef,
+ writer=>'_nhActual');
+has 'converged'=>(is=>'ro', isa=>Num, init_arg=>undef, writer=>'_converged');
+has 'Gpp'=>(is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'Value of projected Greens function');
+has 'waveOperator' => (is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'Wave operator');
+has 'epsilon' => (is=>'lazy', isa=>PDLComplex, init_arg=>undef,
+ documentation=>'Projected dielectric function');
+
+sub _build_Gpp {
+ my $self=shift;
+ my $h = $self->haydock;
+ $h->run unless $h->iteration;
+ my $epsR=$h->epsilonR;
+ my $as=$h->as;
+ my $bcs=$h->bcs;
+ my $min= min($self->nh, $h->iteration);
+ # b0+a1/b1+a2/...
+ # lo debo convertir a
+ # 1-a_0-g0g1b1^2/1-a1-g1g2b2^2/...
+ # entonces bn->1-an y an->-g{n-1}gnbn^2 o -bc_n
+ my ($fn, $n)=lentzCF(1-$as, -$bcs, $min, $self->smallE);
+ #If there are less available coefficients than $self->nh and all
+ #of them were used, there is no remaining work to do, so, converged
+ $self->_converged($n<$min || $h->iteration<=$self->nh);
+ $self->_nhActual($n);
+ my $g0b02=$h->gs->slice("(0)")*$h->b2s->slice("(0)");
+ return $g0b02/($epsR*$fn);
+}
+
+sub _build_waveOperator {
+ my $self=shift;
+ my $greenP=$self->Gpp;
+ my $wave=1/$greenP; #only works along principal directions!!
+ return $wave;
+}
+
+sub _build_epsilon {
+ my $self=shift;
+ my $wave=$self->waveOperator;
+ my $q=$self->haydock->metric->wavenumber;
+ my $q2=$q*$q;
+ my $k=$self->haydock->metric->wavevector;
+ my $k2=($k*$k)->sumover; #inner. my $k2=$k->inner($k); only works on real
+ my $p=$self->haydock->normalizedPolarization;
+ #Note $p->inner($p) might be complex, so is not necessarily 1.
+ my $p2=($p*$p)->sumover;
+ my $pk=($p*$k)->sumover;
+ my $proj=$p2*$k2/$q2 - $pk*$pk/$q2;
+ my $eps=$wave+$proj;
+ return $eps;
+}
+
+__PACKAGE__->meta->make_immutable;
+
+1;
diff --git a/t/00-load.t b/t/00-load.t
index 431a1996..287da97e 100644
--- a/t/00-load.t
+++ b/t/00-load.t
@@ -82,12 +82,12 @@ Photonic::WE::S::Green
Photonic::WE::S::GreenP
Photonic::WE::S::Haydock
Photonic::WE::S::Metric
+Photonic::WE::ST::Field
+Photonic::WE::ST::Green
+Photonic::WE::ST::GreenP
Photonic::WE::ST::Haydock
Photonic::WE::ST::Metric
);
-#Photonic::WE::ST::Field
-#Photonic::WE::ST::Green
-#Photonic::WE::ST::GreenP
foreach(@mods){
diff --git a/t/epst-west.t b/t/epst-west.t
new file mode 100644
index 00000000..5d85f5b4
--- /dev/null
+++ b/t/epst-west.t
@@ -0,0 +1,137 @@
+=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 PDL::NiceSlice;
+use Photonic::Geometry::FromEpsilonTensor;
+use Photonic::WE::ST::Metric;
+use Photonic::WE::ST::Haydock;
+use Photonic::WE::ST::GreenP;
+use Photonic::WE::ST::Green;
+
+use Test::More;
+use lib 't/lib';
+use TestUtils;
+
+#Check epsilontensor for simple 1D system
+#Non-retarded limit
+my ($ea, $eb)=(1+2*i,3+4*i);
+my $f=6/11;
+my $eps=($ea*(zeroes(11,1)->xvals<5)+ $eb*(zeroes(11,1)->xvals>=5))->(*1,*1)*identity(2);
+my $g=Photonic::Geometry::FromEpsilonTensor
+ ->new(epsilon=>$eps);
+my $m=Photonic::WE::ST::Metric->new(
+ geometry=>$g, epsilon=>pdl(1), wavenumber=>pdl(2e-5),
+ wavevector=>pdl([1,0])*2.1e-5);
+my $et=Photonic::WE::ST::Green->new(nh=>10, metric=>$m);
+my $etv=$et->epsilonTensor;
+ok(Cagree($etv->((0),(0)), 1/((1-$f)/$ea+$f/$eb)),
+ "Long. perp. non retarded");
+ok(Cagree($etv->((1),(1)), (1-$f)*$ea+$f*$eb),
+ "Trans. parallel non retarded");
+ok(Cagree($etv->((0),(1)), 0), "xy k perp");
+ok(Cagree($etv->((1),(0)), 0), "yx k perp");
+
+$m=Photonic::WE::ST::Metric->new(
+ geometry=>$g, epsilon=>pdl(1), wavenumber=>pdl(2e-5),
+ wavevector=>pdl([0,1])*2.1e-5);
+$et=Photonic::WE::ST::Green->new(nh=>10, metric=>$m);
+$etv=$et->epsilonTensor;
+ok(Cagree($etv->((0),(0)), 1/((1-$f)/$ea+$f/$eb)),
+ "Trans. perp. non retarded");
+ok(Cagree($etv->((1),(1)), (1-$f)*$ea+$f*$eb),
+ "Long. parallel non retarded");
+ok(Cagree($etv->((0),(1)), 0), "xy k parallel");
+ok(Cagree($etv->((1),(0)), 0), "yx k parallel");
+
+#Compare to epsilon from transfer matrix.
+#Construct normal incidence transfer matrix
+($ea, $eb)=(r2C(1),r2C(2));
+$eps=r2C($ea*(zeroes(11,1)->xvals<5)+ $eb*(zeroes(11)->xvals>=5))->(*1,*1)*identity(2);
+$g=Photonic::Geometry::FromEpsilonTensor
+ ->new(epsilon=>$eps, L=>pdl(1,1));
+my ($na, $nb)=map {sqrt($_)} ($ea, $eb);
+my $q=1.2;
+my ($ka,$kb)=map {$q*$_} ((1-$f)*$na, $f*$nb); #Multiply by length also
+my $ma=pdl([cos($ka), -sin($ka)/$na],[$na*sin($ka), cos($ka)]);
+my $mb=pdl([cos($kb), -sin($kb)/$nb],[$nb*sin($kb), cos($kb)]);
+my $mt=($ma->(:,*1)*$mb->transpose->(:,:,*1))->sumover;
+#Solve exact dispersion relation
+my $cospd=($mt->((0),(0))+$mt->((1),(1)))/2;
+my $sinpd=sqrt(1-$cospd**2);
+my $pd=log($cospd+i()*$sinpd)/i;
+warn "Bloch vector not real, $pd" unless $pd->im->abs < 1e-7;
+$pd=$pd->re;
+#epsilon from transfer matrix
+my $epstm=($pd/$q)**2;
+#epsilon from photonic
+$m=Photonic::WE::ST::Metric->new(
+ geometry=>$g, epsilon=>pdl(1), wavenumber=>pdl($q),
+ wavevector=>pdl([$pd,0]));
+$et=Photonic::WE::ST::Green->new(nh=>1000, metric=>$m,
+ reorthogonalize=>1);
+$etv=$et->epsilonTensor->((1),(1));
+ok(Cagree($epstm, $etv), "Epsilon agrees with transfer matrix");
+
+#Compare to epsilon from transfer matrix with complex metric.
+#Construct normal incidence transfer matrix
+($ea, $eb)=(r2C(1),2+1*i);
+$eps=($ea*(zeroes(11,1)->xvals<5)+ $eb*(zeroes(11)->xvals>=5))->(*1,*1)*identity(2);
+$g=Photonic::Geometry::FromEpsilonTensor->new(epsilon=>$eps, L=>pdl(1,1));
+($na, $nb)=map {sqrt($_)} ($ea, $eb);
+$q=1.2;
+($ka,$kb)=map {$q*$_} ((1-$f)*$na, $f*$nb); #Multiply by length also
+$ma=pdl([cos($ka), -sin($ka)/$na],[$na*sin($ka), cos($ka)]);
+$mb=pdl([cos($kb), -sin($kb)/$nb],[$nb*sin($kb), cos($kb)]);
+$mt=($ma->(:,*1)*$mb->transpose->(:,:,*1))->sumover;
+#Solve exact dispersion relation
+$cospd=($mt->((0),(0))+$mt->((1),(1)))/2;
+$sinpd=sqrt(1-$cospd**2);
+$pd=log($cospd+i()*$sinpd)/i;
+#epsilon from transfer matrix
+$epstm=($pd/$q)**2;
+#epsilon from photonic
+$m=Photonic::WE::ST::Metric->new(
+ geometry=>$g, epsilon=>pdl(1), wavenumber=>pdl($q),
+ wavevector=>pdl([$pd,0]));
+$et=Photonic::WE::ST::Green->new(nh=>1000, metric=>$m,
+ reorthogonalize=>1);
+$etv=$et->epsilonTensor->((1),(1));
+ok(Cagree($epstm, $etv, 1e-4), "Epsilon agrees with transfer matrix. Complex case.");
+
+my $h=Photonic::WE::ST::Haydock->new(nh=>1000, metric=>$m,
+ polarization=>r2C(pdl [0,1]), reorthogonalize=>1);
+$et=Photonic::WE::ST::GreenP->new(nh=>1000, haydock=>$h);
+$etv=$et->epsilon;
+ok(Cagree($epstm, $etv, 1e-4), "Projected eps agrees with trans mat. Complex case.");
+
+done_testing();
diff --git a/t/field-west.t b/t/field-west.t
new file mode 100644
index 00000000..8f5c2f5e
--- /dev/null
+++ b/t/field-west.t
@@ -0,0 +1,75 @@
+=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::WE::ST::Haydock;
+use Photonic::WE::ST::Metric;
+use Photonic::WE::ST::Field;
+
+use Test::More tests => 2;
+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)->slice("*1,*1"); # *identity(1);
+my $gl=Photonic::Geometry::FromB->new(B=>$B); #long
+my $ml=Photonic::WE::ST::Metric->new(geometry=>$gl, epsilon=>pdl(1),
+ wavenumber=>pdl(1), wavevector=>pdl([0.01]));
+my $haydock=Photonic::WE::ST::Haydock->new(
+ metric=>$ml, nh=>10, keepStates=>1, polarization=>pdl([1])->r2C,
+ epsilon=>$epsilon);
+my $flo=Photonic::WE::ST::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)->slice("*1,*1")*identity(2);
+my $gt=Photonic::Geometry::FromB->new(B=>$Bt); #trans
+my $mt=Photonic::WE::ST::Metric->new(geometry=>$gt, epsilon=>pdl(1),
+ wavenumber=>pdl(0.001), wavevector=>pdl([0,0.0001]));
+my $nt=Photonic::WE::ST::Haydock->new(
+ metric=>$mt, nh=>10, keepStates=>1, polarization=>pdl([1,0])->r2C,
+ epsilon=>$epsilont);
+my $fto=Photonic::WE::ST::Field->new(haydock=>$nt, nh=>10);
+my $ftv=$fto->field;
+my $ftx=r2C(pdl [1, 0]);
+ok(Cagree($ftv, $ftx), "1D trans field");
diff --git a/t/oneh-west.t.t b/t/oneh-west.t
similarity index 100%
rename from t/oneh-west.t.t
rename to t/oneh-west.t