diff --git a/HW3/homework3/src/compute_temperature.cc b/HW3/homework3/src/compute_temperature.cc
index 4932cbc..63e514b 100644
--- a/HW3/homework3/src/compute_temperature.cc
+++ b/HW3/homework3/src/compute_temperature.cc
@@ -1,72 +1,67 @@
 #include "compute_temperature.hh"
 #include "fft.hh"
 #include "material_point.hh"
 #include <cmath>
 
 /* -------------------------------------------------------------------------- */
 
-ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {}
+ComputeTemperature::ComputeTemperature(Real dt,double boundary) : dt(dt),boundary(boundary) {}
 
 void ComputeTemperature::compute(System& system)
 {
   UInt npnt=system.getNbParticles();
   UInt size=sqrt(npnt);
-  Real L=2.0;
+  Real L=1.0;
   Real k=1;
   Matrix<complex> temp(size),heatf(size),dtemp(size);
   Matrix<complex> Ftemp(size),Fheatf(size),Fdtemp(size);
   Matrix<std::complex<int>> freq=FFT::computeFrequencies(size);
 
 
   for (UInt i = 0; i < npnt; i++)
   {
     //MaterialPoint &par=system.getParticle(i);
     MaterialPoint &par = static_cast<MaterialPoint &>(system.getParticle(i));
     temp(i)=par.getTemperature();
     heatf(i)=par.getHeatRate();
   }
   auto *vtemp=temp.data();
 
 
   Ftemp=FFT::transform(temp);
   Fheatf=FFT::transform(heatf);
   Fdtemp=Fheatf-k*4*M_PI*M_PI/(L*L)*freq*Ftemp;
   dtemp=FFT::itransform(Fdtemp);
   temp=temp+dt*dtemp;
-
+  UInt N=size;
   for (UInt i = 0; i < npnt; i++)
   {
     //MaterialPoint &par=system.getParticle(i);
     MaterialPoint &par = static_cast<MaterialPoint &>(system.getParticle(i));
-    if(i/size==0||i%size==0||i%size==size-1||i/size==size-1)
+    if(!std::isnan(boundary))
     {
-    temp(i)=0;
+      if((i%N)==0 || (i%N)==(N-1) || (i/N)==0 || (i/N)==(N-1))
+      {
+        temp(i).real(boundary);
+      }
+    }
+
     par.getTemperature()=temp(i).real();
-   } else
-   par.getTemperature()=temp(i).real();
-   {
   }
-}
-}
-
-
-    //
-    //;
-
 
 
   /*auto *vtemp=temp.data(),*vheatf=heatf.data(),*vFtemp=Ftemp.data(),
   *vFdtemp=Fdtemp.data(),*vFheatf=Fheatf.data();
   auto *vfreq=freq.data();
   //test=Fdtemp*Fheatf;
 
   for(int i=0;i<npnt;i++)
   {
     *(vFdtemp+i)=*(vFheatf)-k*4*M_PI*M_PI/npnt*(*(vfreq+i))*(*(vFtemp+i));
   }*/
   //Fdtemp=Fheatf-k*4*M_PI*M_PI/npnt*FFT::computeFrequencies(sqrt(npnt))*Ftemp;
 
 
-
+}
 
 /* -------------------------------------------------------------------------- */
diff --git a/HW3/homework3/src/compute_temperature.hh b/HW3/homework3/src/compute_temperature.hh
index a61a43d..72a8cc8 100644
--- a/HW3/homework3/src/compute_temperature.hh
+++ b/HW3/homework3/src/compute_temperature.hh
@@ -1,22 +1,24 @@
 #ifndef __COMPUTE_TEMPERATURE__HH__
 #define __COMPUTE_TEMPERATURE__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "compute.hh"
+#include <limits>
 
 //! Compute contact interaction between ping-pong balls
 class ComputeTemperature : public Compute {
 
   // Virtual implementation
 public:
-    ComputeTemperature(Real dt);
+    ComputeTemperature(Real dt,double bboundary=std::numeric_limits<double>::quiet_NaN());
   //! Penalty contact implementation
   void compute(System& system) override;
 
 private:
     Real dt;
+    double boundary;
 
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__COMPUTE_TEMPERATURE__HH__
diff --git a/HW3/homework3/src/fft.hh b/HW3/homework3/src/fft.hh
index 2fbb901..e490592 100644
--- a/HW3/homework3/src/fft.hh
+++ b/HW3/homework3/src/fft.hh
@@ -1,110 +1,137 @@
 #ifndef FFT_HH
 #define FFT_HH
 /* ------------------------------------------------------ */
 #include "matrix.hh"
 #include "my_types.hh"
 #include <fftw3.h>
 /* ------------------------------------------------------ */
 
 struct FFT {
 
   static Matrix<complex> transform(Matrix<complex>& m);
   static Matrix<complex> itransform(Matrix<complex>& m);
 
   static Matrix<std::complex<int>> computeFrequencies(int size);
 
 };
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::transform(Matrix<complex>& m_in)
 {
   fftw_complex *in, *out;
   fftw_plan p;
   UInt size=m_in.size();
   auto *val_in=m_in.data();
 
   in = new fftw_complex[size*size];
   out = new fftw_complex[size*size];
 
   p = fftw_plan_dft_2d(size,size, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
 
   for(int i=0;i<size*size;i++)
   {
     in[i][0]=(val_in+i)->real();
     in[i][1]=(val_in+i)->imag();
   }
 
   fftw_execute(p);
 
   Matrix<complex> m_out(size);
 
   auto *val_out=m_out.data();
 
   for(int i=0;i<size*size;i++)
   {
     *(val_out+i)=std::complex<double>(out[i][0],out[i][1]);
   }
 
   fftw_destroy_plan(p);
   delete [] in;
   delete [] out;
 
   return(m_out);
 }
 
 /* ------------------------------------------------------ */
 
 inline Matrix<complex> FFT::itransform(Matrix<complex>& m_in)
 {
   fftw_complex *in, *out;
   fftw_plan p;
   UInt size=m_in.size();
   auto *val_in=m_in.data();
 
   in = new fftw_complex[size*size];
   out = new fftw_complex[size*size];
 
   p = fftw_plan_dft_2d(size,size, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
 
   for(int i=0;i<size*size;i++)
   {
     in[i][0]=(val_in+i)->real();
     in[i][1]=(val_in+i)->imag();
   }
 
   fftw_execute(p);
 
   Matrix<complex> m_out(size);
 
   auto *val_out=m_out.data();
 
   for(int i=0;i<size*size;i++)
   {
     *(val_out+i)=std::complex<double>(out[i][0],out[i][1])/(double)(size*size);
   }
 
   fftw_destroy_plan(p);
   delete [] in;
   delete [] out;
 
   return(m_out);
 }
 
 /* ------------------------------------------------------ */
 
 
 /* ------------------------------------------------------ */
 
 inline Matrix<std::complex<int>> FFT::computeFrequencies(int size)
 {
   Matrix<std::complex<int>> freq(size);
   auto *value=freq.data();
   for(int i=0;i<size*size;i++)
   {
-    *(value+i)=std::complex<int>((i%size)*(i%size)+(i/size)*(i/size),0);
+    if(i%size<(size/2+size%2))
+    {
+      (value+i)->real((i%size)*(i%size));
+    }
+    else
+    {
+      (value+i)->real((i%size-size)*(i%size-size));
+    }
+  }
+  for(int i=0;i<size*size;i++)
+  {
+    if(i/size<size/2+size%2)
+    {
+      (value+i)->real((value+i)->real()+(i/size)*(i/size));
+    }
+    else
+    {
+      (value+i)->real((value+i)->real()+(i/size-size)*(i/size-size));
+    }
+  }
+  for(int i=0;i<size*size;i++)
+  {
+    (value+i)->imag((value+i)->real());
   }
+
+  /*for(int i=0;i<size*size;i++)
+  {
+    *(value+i)=std::complex<int>((i%size)*(i%size)+(i/size)*(i/size),0);
+  }*/
   return(freq);
 }
 
 #endif  // FFT_HH
diff --git a/HW3/homework3/src/main.cc b/HW3/homework3/src/main.cc
index a50a259..97cc471 100644
--- a/HW3/homework3/src/main.cc
+++ b/HW3/homework3/src/main.cc
@@ -1,62 +1,71 @@
 #include "compute_gravity.hh"
 #include "compute_verlet_integration.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "my_types.hh"
 #include "ping_pong_balls_factory.hh"
 #include "material_points_factory.hh"
 #include "planets_factory.hh"
 #include "system.hh"
 /* -------------------------------------------------------------------------- */
 #include <cstdlib>
 #include <iostream>
 #include <sstream>
 /* -------------------------------------------------------------------------- */
 
 int main(int argc, char** argv) {
-  if (argc != 6) {
+  if (argc < 6) {
     std::cout << "Usage: " << argv[0]
               << " nsteps dump_freq input.csv particle_type timestep"
               << std::endl;
     std::cout << "\tparticle type can be: planet, ping_pong, material_point" << std::endl;
     std::exit(EXIT_FAILURE);
   }
 
 
   // the number of steps to perform
   Real nsteps;
   std::stringstream(argv[1]) >> nsteps;
   // freq to dump
   int freq;
   std::stringstream(argv[2]) >> freq;
   // init file
   std::string filename = argv[3];
   // particle type
   std::string type = argv[4];
   // timestep
   Real timestep;
   std::stringstream(argv[5]) >> timestep;
 
   if (type == "planet")
     PlanetsFactory::getInstance();
   else if (type == "ping_pong")
     PingPongBallsFactory::getInstance();
   else if (type == "material_point")
-    MaterialPointsFactory::getInstance();  
+  {
+    MaterialPointsFactory::getInstance();
+    if(argc==7)
+    {
+      MaterialPointsFactory &factory = static_cast<MaterialPointsFactory &>(ParticlesFactoryInterface::getInstance());
+      double boundary;
+      std::stringstream(argv[6])>>boundary;
+      factory.setboundary(boundary);
+    }
+  }
   else {
     std::cout << "Unknown particle type: " << type << std::endl;
     std::exit(EXIT_FAILURE);
   }
 
+
   ParticlesFactoryInterface& factory = ParticlesFactoryInterface::getInstance();
 
   SystemEvolution& evol = factory.createSimulation(filename, timestep);
-
   evol.setNSteps(nsteps);
   evol.setDumpFreq(freq);
 
   evol.evolve();
 
 
   return EXIT_SUCCESS;
 }
diff --git a/HW3/homework3/src/material_points_factory.cc b/HW3/homework3/src/material_points_factory.cc
index f44e30d..bfc13ae 100644
--- a/HW3/homework3/src/material_points_factory.cc
+++ b/HW3/homework3/src/material_points_factory.cc
@@ -1,47 +1,56 @@
 #include "material_points_factory.hh"
 #include "compute_temperature.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "material_point.hh"
 #include <cmath>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 std::unique_ptr<Particle> MaterialPointsFactory::createParticle() {
   return std::make_unique<MaterialPoint>();
 }
 
 /* -------------------------------------------------------------------------- */
 
 SystemEvolution&
 MaterialPointsFactory::createSimulation(const std::string& fname,
                                         Real timestep) {
 
+
   this->system_evolution =
       std::make_unique<SystemEvolution>(std::make_unique<System>());
 
   CsvReader reader(fname);
   reader.read(this->system_evolution->getSystem());
 
   // check if it is a square number
   auto N = this->system_evolution->getSystem().getNbParticles();
   int side = std::sqrt(N);
   if (side * side != N)
     throw std::runtime_error("number of particles is not square");
 
-  auto temperature = std::make_shared<ComputeTemperature>(timestep);
-  this->system_evolution->addCompute(temperature);
+  if(!std::isnan(boundary))
+  {
+    auto temperature = std::make_shared<ComputeTemperature>(timestep,boundary);
+    this->system_evolution->addCompute(temperature);
+  }
+  else
+  {
+    auto temperature = std::make_shared<ComputeTemperature>(timestep);
+    this->system_evolution->addCompute(temperature);
+  }
 
   return *system_evolution;
 }
 
 /* -------------------------------------------------------------------------- */
 
 ParticlesFactoryInterface& MaterialPointsFactory::getInstance() {
   if (not ParticlesFactoryInterface::factory)
     ParticlesFactoryInterface::factory = new MaterialPointsFactory;
 
   return *factory;
 }
 
 /* -------------------------------------------------------------------------- */
diff --git a/HW3/homework3/src/material_points_factory.hh b/HW3/homework3/src/material_points_factory.hh
index 85793cd..9ced37c 100644
--- a/HW3/homework3/src/material_points_factory.hh
+++ b/HW3/homework3/src/material_points_factory.hh
@@ -1,30 +1,39 @@
 #ifndef __MATERIAL_POINTS_FACTORY__HH__
 #define __MATERIAL_POINTS_FACTORY__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "particles_factory_interface.hh"
 #include "material_point.hh"
+#include <limits>
 /* -------------------------------------------------------------------------- */
 
 //! Factory for material points
 class MaterialPointsFactory : public ParticlesFactoryInterface {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 private:
   MaterialPointsFactory() = default;
 
+  double boundary=std::numeric_limits<double>::quiet_NaN();
+
+
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
+
   SystemEvolution& createSimulation(const std::string& fname,
 				    Real timestep) override;
 
   std::unique_ptr<Particle> createParticle() override;
 
   static ParticlesFactoryInterface& getInstance();
+
+  void setboundary(double boundary){this->boundary=boundary;}
+
+  void getboundary(){std::cout<<boundary;}
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__MATERIAL_POINTS_FACTORY__HH__
diff --git a/HW3/homework3/src/particles_factory_interface.hh b/HW3/homework3/src/particles_factory_interface.hh
index b04c38b..8190be6 100644
--- a/HW3/homework3/src/particles_factory_interface.hh
+++ b/HW3/homework3/src/particles_factory_interface.hh
@@ -1,39 +1,41 @@
 #ifndef __PARTICLES_FACTORY_INTERFACE__HH__
 #define __PARTICLES_FACTORY_INTERFACE__HH__
 
 /* -------------------------------------------------------------------------- */
 #include "system_evolution.hh"
 /* -------------------------------------------------------------------------- */
 
 //! Abstract factory defining interface
 class ParticlesFactoryInterface {
   // Constructors/Destructors
 protected:
   //! Instance constructor (protected)
   ParticlesFactoryInterface() = default;
 
 public:
   virtual ~ParticlesFactoryInterface() = default;
 
   // Methods
 public:
   //! Create a whole simulation from file
   virtual SystemEvolution& createSimulation(const std::string& fname,
                                             Real timestep) = 0;
   //! Create a new particle
   virtual std::unique_ptr<Particle> createParticle() = 0;
 
   //! Get singleton instance
   static ParticlesFactoryInterface& getInstance();
 
+
+
   // Members
 protected:
   std::vector<Particle*> list_particles;
   std::unique_ptr<SystemEvolution> system_evolution = nullptr;
 
   // Standard pointer because constructor is protected (cannot use make_unique)
   static ParticlesFactoryInterface* factory;
 };
 
 /* -------------------------------------------------------------------------- */
 #endif  //__PARTICLES_FACTORY_INTERFACE__HH__
diff --git a/HW3/homework3/src/test_computetemperature.cc b/HW3/homework3/src/test_computetemperature.cc
index a74dc1d..e725c65 100644
--- a/HW3/homework3/src/test_computetemperature.cc
+++ b/HW3/homework3/src/test_computetemperature.cc
@@ -1,100 +1,100 @@
 #include "my_types.hh"
 #include "fft.hh"
 #include "system.hh"
 #include "material_points_factory.hh"
 #include <gtest/gtest.h>
 #include "compute_temperature.hh"
-#include <fstream>
 
 /*****************************************************************/
 TEST(CT, UniformT) {
-  UInt N = 512;
+  UInt N = 64;
 
-  Real nsteps=20,initT=300,timestep=0.000001;
+  Real nsteps=20,initT=300,timestep=0.0001;
   System system;
   double dx=2.0/(N-1);
 
   MaterialPointsFactory::getInstance();
 
   for(int i=0;i<N*N;i++)
   {
     std::stringstream sstr;
     sstr<<(i%N)*dx-1<<" "<<(i/N)*dx-1<<" "<<0<<" ";
     sstr<<0<<" "<<0<<" "<<0<<" ";
     sstr<<0<<" "<<0<<" "<<0<<" ";
     sstr<<0<<" "<<initT<<" "<<0<<" ";
     auto p = ParticlesFactoryInterface::getInstance().createParticle();
     sstr >> *p;
     system.addParticle(std::move(p));
   }
 
-  ComputeTemperature temperature(timestep);
+  ComputeTemperature temperature(0.001);
 
     for (UInt i = 0; i < nsteps; ++i)
     {
         temperature.compute(system);
     }
 
     for(int i=0;i<N*N;i++)
     {
         MaterialPoint &p = static_cast<MaterialPoint &>(system.getParticle(i));
         //std::cout<<p.getTemperature()<<std::endl;
         ASSERT_NEAR(p.getTemperature(),initT,1e-10);
         ASSERT_NEAR(p.getHeatRate(),0,1e-10);
     }
 }
 
 /*****************************************************************/
+TEST(CT, NonUniformT)
+{
+  UInt N = 16;
 
+  Real nsteps=400,T1=20,T2=25,initT,timestep=0.001;
+  System system;
+  double dx=2.0/(N-1);
 
-  TEST(CT, Heatflux) {
-    UInt N = 128;
-
-    Real nsteps=20,initT=0,timestep=0.000001, initheatf;
-    Real F1=1, F2=-1;
-    System system;
-    double dx=2.0/(N-1);
-
-    MaterialPointsFactory::getInstance();
+  MaterialPointsFactory::getInstance();
 
-    for(int i=0;i<N*N;i++)
+  for(int i=0;i<N*N;i++)
+  {
+    if(((i%N)*dx-1)<=0)
     {
-      if(i==N/4)
-      {
-        initheatf=-N;
-      }
-      else if (i==3*N/4)
-      {
-        initheatf=-N;
-      }
-      else
-      {
-       initheatf=0;
-      }
-      std::stringstream sstr;
-      sstr<<(i%N)*dx-1<<" "<<(i/N)*dx-1<<" "<<0<<" ";
-      sstr<<0<<" "<<0<<" "<<0<<" ";
-      sstr<<0<<" "<<0<<" "<<0<<" ";
-      sstr<<0<<" "<<initT<<" "<<initheatf<<" ";
-      auto p = ParticlesFactoryInterface::getInstance().createParticle();
-      sstr >> *p;
-      system.addParticle(std::move(p));
+      initT=T1;
     }
-
-    ComputeTemperature temperature(timestep);
-
-    for (UInt i = 0; i < nsteps; ++i)
+    else
     {
-        temperature.compute(system);
+      initT=T2;
     }
+    std::stringstream sstr;
+    sstr<<(i%N)*dx-1<<" "<<(i/N)*dx-1<<" "<<0<<" ";
+    sstr<<0<<" "<<0<<" "<<0<<" ";
+    sstr<<0<<" "<<0<<" "<<0<<" ";
+    sstr<<0<<" "<<initT<<" "<<0<<" ";
+    auto p = ParticlesFactoryInterface::getInstance().createParticle();
+    sstr >> *p;
+    system.addParticle(std::move(p));
+  }
 
-    std::ifstream fin;
-    fin.open("ouput.csv");
-    double temp;
-     for (int i=0;i<N*N;i++)
-    {
-        fin>>temp;
-        MaterialPoint &p = static_cast<MaterialPoint &>(system.getParticle(i));
-        ASSERT_NEAR(p.getTemperature(),temp,1e-4);
-    }
+  ComputeTemperature temperature(timestep);
+
+  for (UInt i = 0; i < nsteps; ++i)
+  {
+      temperature.compute(system);
   }
+
+  for(int i=0;i<N*N;i++)
+  {
+      MaterialPoint &p = static_cast<MaterialPoint &>(system.getParticle(i));
+      //std::cout<<p.getTemperature()<<std::endl;
+      ASSERT_NEAR(p.getTemperature(),(T1+T2)/2,1);
+  }
+}
+
+/*****************************************************************/
+/*TEST(CT, NonUniformT)
+{
+  UInt N = 128;
+
+  Real nsteps=20,T1=300,T2=200,initT,timestep=0.000001;
+  System system;
+  double dx=2.0/(N-1);
+}*/