Initial commitHEADmaster

1 files changed, 228 insertions, 0 deletions

diff --git a/fresnel/propagate.C b/fresnel/propagate.C
new file mode 100644
index 0000000..0dae533
--- /dev/null
+++ b/fresnel/propagate.C
@@ -0,0 +1,228 @@
+/*
+   This file is part of LPIC++, a particle-in-cell code for
+   simulating the interaction of laser light with plasma.
+
+   Copyright (C) 1994-1997 Roland Lichters
+
+   LPIC++ 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 2
+   of the License, 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+*/
+
+#include <propagate.h>
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+propagate::propagate(parameter &p)
+    : zeit(p)
+{
+  sprintf( errname, "%s/error", p.path );
+  static error_handler bob("propagate::Constructor",errname);
+
+  dt    = 1.0 / p.prop.spp;
+  dw    = 1.0 / ( p.prop.time_stop - p.prop.time_start ) / 50;
+  wmax  = 5;
+  t     = p.prop.time_start;
+  T     = p.pulse.duration;
+
+  reflex_avg = 0;
+
+  strcpy(p.errname,errname);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void propagate::loop(parameter &p)
+{
+  static error_handler bob("propagate::loop",errname);
+
+  out_int(p);      // plot the integrand
+
+  for( t = p.prop.time_start, t_avg=0; t <= p.prop.time_stop+dt; t += dt, t_avg++ )
+    {
+      w = 0.5*dw; reflex = 0;
+      do
+	{
+	  r = reflectivity( w, p.plasma.wp, p.pulse.angle_rad, p.pulse.polarization );
+	  phi = phase( w, p.plasma.wp, p.pulse.angle_rad, p.pulse.polarization );
+	  integrand = r * ( dif((w-1)*PI*T) - dif((w+1)*PI*T) ) * sin( w*PI*(2*t-T)-phi );
+	  reflex += dw * T * integrand;
+
+	  w += dw;
+	}
+      while( w<wmax );
+
+      out(p);    // plot the flux (f(t)/f0)^2
+
+      reflex_avg += 2 * pow( reflex, 2 ) / p.prop.spp;
+
+                 // plot the ratio of time averaged fluxes once per period
+      if ( t_avg==p.prop.spp ) { out_avg(p); t_avg=0; reflex_avg=0; }
+
+      zeit.proc();              // update clock
+    }
+
+  zeit.exit();
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+double propagate::reflectivity( double w, double wp, double angle, int polarization )
+{
+  static error_handler bob("propagate::reflectivity",errname);
+
+  if (w*cos(angle)>wp) {
+
+    double n = sqrt( 1 - pow(wp/w,2) );
+
+    rs = cos(angle) - sqrt( n*n-sin(angle)*sin(angle) );
+    rs /= ( cos(angle) + sqrt( n*n-sin(angle)*sin(angle) ) );
+
+    rp = sqrt( 1-pow(sin(angle)/n,2) ) - n*cos(angle);
+    rp /= ( sqrt( 1-pow(sin(angle)/n,2) ) + n*cos(angle) );
+  }
+  else {
+    rs = rp = 1;
+  }
+
+  if (polarization==1) return rs;
+  else                 return rp;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+double propagate::phase( double w, double wp, double angle, int polarization )
+{
+  static error_handler bob("propagate::phase",errname);
+  double brev = w*cos(angle)/wp;
+
+
+  if (brev>=1) {
+
+    double n = sqrt( 1 - pow(wp/w,2) );
+
+    phis = phip = 0;
+  }
+  else {
+
+    // from Born and Wolf, page 49
+    phis = - 2 * atan( sqrt( 1-brev*brev ) / brev );
+
+    // from Born and Wolf, page 49
+    phip = + 2 * atan( wp/w * cos(angle) * ( pow(w/wp,2)-1 ) / sqrt( 1-brev*brev ) );
+  }
+
+  if (polarization==1) return phis;
+  else                 return phip;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+double propagate::dif( double x )
+{
+  static error_handler bob("propagate::dif",errname);
+
+  if ( x==0 ) return 1;
+  else return sin(x)/x;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void propagate::out( parameter &p )
+{
+  static error_handler bob("propagate::out",errname);
+
+  char filename[filename_size];
+
+  sprintf(filename,"%s/reflex", p.path );
+
+  reflex_file.open(filename,ios::app);
+  if (!reflex_file) bob.error( "cannot open file", filename );
+
+  reflex_file.precision( 3 );
+  reflex_file.setf( ios::showpoint | ios::scientific );
+
+  reflex_file << setw(10)  << t  << setw(10) << pow( p.pulse.amplitude*reflex,2 ) << endl;
+
+  reflex_file.close();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void propagate::out_avg( parameter &p )
+{
+  static error_handler bob("propagate::out_avg",errname);
+
+  char filename[filename_size];
+
+  sprintf(filename,"%s/reflex-avg", p.path );
+
+  reflex_file_avg.open(filename,ios::app);
+  if (!reflex_file_avg) bob.error( "cannot open file", filename );
+
+  reflex_file_avg.precision( 3 );
+  reflex_file_avg.setf( ios::showpoint | ios::scientific );
+
+  reflex_file_avg << setw(10)  << t  << setw(10) << reflex_avg << endl;
+
+  reflex_file_avg.close();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void propagate::out_int( parameter &p )
+{
+  static error_handler bob("propagate::out_int",errname);
+
+  char filename[filename_size];
+
+  sprintf(filename,"%s/reflex-int", p.path );
+
+  reflex_file_int.open(filename,ios::app);
+  if (!reflex_file_int) bob.error( "cannot open file", filename );
+
+  reflex_file_int.precision( 3 );
+  reflex_file_int.setf( ios::showpoint | ios::scientific );
+
+  reflex_file_int << "#"
+                  << setw(12)  << "w"     << setw(13) << "r(w)"
+		  << setw(13)  << "phase"
+                  << setw(13)  << "dif()" << setw(13) << "r(w)*dif()" << endl;
+
+  w = 0; reflex = 0;
+  do
+    {
+      r = reflectivity( w, p.plasma.wp, p.pulse.angle_rad, p.pulse.polarization );
+      phi = phase( w, p.plasma.wp, p.pulse.angle_rad, p.pulse.polarization );
+      integrand = r * ( dif( (w-1)*PI*T ) - dif( (w+1)*PI*T) );
+
+      reflex_file_int << setw(13) << w
+                      << setw(13) << r
+		      << setw(13) << phi
+		      << setw(13) << dif( (w-1)*PI*T ) - dif( (w+1)*PI*T)
+		      << setw(13) << integrand << endl;
+
+      w += 5*dw;
+    }
+  while( w<2 );
+
+  reflex_file_int.close();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+//EOF
+
+
+