diff options
Diffstat (limited to 'fresnel/propagate.C')
| -rw-r--r-- | fresnel/propagate.C | 228 |
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 + + + |