summaryrefslogtreecommitdiff
path: root/post/src/spacetime.C
diff options
context:
space:
mode:
Diffstat (limited to 'post/src/spacetime.C')
-rw-r--r--post/src/spacetime.C1087
1 files changed, 1087 insertions, 0 deletions
diff --git a/post/src/spacetime.C b/post/src/spacetime.C
new file mode 100644
index 0000000..1f22b3b
--- /dev/null
+++ b/post/src/spacetime.C
@@ -0,0 +1,1087 @@
+/*
+ This file is part of LPIC++, a particle-in-cell code for
+ simulating the interaction of laser light with plasma.
+
+ Copyright (C) 2002 Andreas Kemp
+ 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.
+*/
+
+//////////////////////////////////////////////////////////////////////////////////////////
+//
+// postprocessor for lpic++
+//
+//////////////////////////////////////////////////////////////////////////////////////////
+// changes by A.Kemp, 2002 denoted by ##
+#include <spacetime.h>
+
+spacetime::spacetime( parameter &p )
+ : input(p),
+ ft( input.periods_x, input.cells_per_wl, 1 ),
+ ft2d( input.Q_kw, input.periods_t, input.steps_per_period, 1,
+ input.periods_x, input.cells_per_wl, 1 )
+{
+ sprintf( errname, "%s/error", p.output_path );
+ static error_handler bob("spacetime::Constructor",errname);
+
+ input_path = new char [ filename_size ];
+ strcpy(input_path,p.file_path);
+ output_path = new char [ filename_size ];
+ strcpy(output_path,p.output_path);
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+input_spacetime::input_spacetime( parameter &p )
+ : rf()
+{
+ strcpy( errname, p.errname );
+ static error_handler bob("input_spacetime::Constructor",errname);
+ char fname[filename_size];
+
+ rf.openinput( p.read_filename );
+
+ t_start = atoi( rf.setget( "&spacetime", "t_start" ) );
+ t_stop = atoi( rf.setget( "&spacetime", "t_stop" ) );
+ x_start = atof( rf.setget( "&spacetime", "x_start" ) );
+ x_stop = atof( rf.setget( "&spacetime", "x_stop" ) );
+ x_offset = atof( rf.setget( "&spacetime", "x_offset" ) );
+
+ periods_x = (int) floor( x_stop - x_start + 0.5 );
+ periods_t = (int) floor( t_stop - t_start + 0.5 );
+
+ average = atoi( rf.setget( "&spacetime", "smooth" ) );
+ size = atoi( rf.setget( "&spacetime", "imagesize" ) );
+ contour_1 = atof( rf.setget( "&spacetime", "contour_1" ) );
+ contour_2 = atof( rf.setget( "&spacetime", "contour_2" ) );
+ contour_3 = atof( rf.setget( "&spacetime", "contour_3" ) );
+
+ Q_kw = atoi( rf.setget( "&spacetime", "Q_kw" ) );
+ Q_kt = atoi( rf.setget( "&spacetime", "Q_kt" ) );
+
+ Q_de = atoi( rf.setget( "&spacetime", "Q_de" ) );
+ Q_di = atoi( rf.setget( "&spacetime", "Q_di" ) );
+ Q_jx = atoi( rf.setget( "&spacetime", "Q_jx" ) );
+ Q_jy = atoi( rf.setget( "&spacetime", "Q_jy" ) );
+ Q_jz = atoi( rf.setget( "&spacetime", "Q_jz" ) );
+ Q_ex = atoi( rf.setget( "&spacetime", "Q_ex" ) );
+ Q_ey = atoi( rf.setget( "&spacetime", "Q_ey" ) );
+ Q_ez = atoi( rf.setget( "&spacetime", "Q_ez" ) );
+ // Q_bx = atoi( rf.setget( "&spacetime", "Q_bx" ) );
+ Q_bx = 0; // ##
+ Q_by = atoi( rf.setget( "&spacetime", "Q_by" ) );
+ Q_bz = atoi( rf.setget( "&spacetime", "Q_bz" ) );
+ Q_edens = atoi( rf.setget( "&spacetime", "Q_edens" ) );
+ Q_de_fi = atoi( rf.setget( "&spacetime", "Q_de_fi" ) );
+ Q_de_ii = atoi( rf.setget( "&spacetime", "Q_de_ii" ) );
+
+ C_kw = atof( rf.setget( "&spacetime", "C_kw" ) );
+ C_kt = atof( rf.setget( "&spacetime", "C_kt" ) );
+
+ K_cut = atof( rf.setget( "&spacetime", "K_cut" ) );
+ W_cut = atof( rf.setget( "&spacetime", "W_cut" ) );
+
+ C_de = atof( rf.setget( "&spacetime", "C_de" ) );
+ C_di = atof( rf.setget( "&spacetime", "C_di" ) );
+ C_jx = atof( rf.setget( "&spacetime", "C_jx" ) );
+ C_jy = atof( rf.setget( "&spacetime", "C_jy" ) );
+ C_jz = atof( rf.setget( "&spacetime", "C_jz" ) );
+ C_ex = atof( rf.setget( "&spacetime", "C_ex" ) );
+ C_ey = atof( rf.setget( "&spacetime", "C_ey" ) );
+ C_ez = atof( rf.setget( "&spacetime", "C_ez" ) );
+ // C_bx = atof( rf.setget( "&spacetime", "C_bx" ) );
+ C_bx = 0; // ##
+ C_by = atof( rf.setget( "&spacetime", "C_by" ) );
+ C_bz = atof( rf.setget( "&spacetime", "C_bz" ) );
+ C_edens = atof( rf.setget( "&spacetime", "C_edens" ) );
+ C_de_fi = atof( rf.setget( "&spacetime", "C_de_fi" ) );
+ C_de_ii = atof( rf.setget( "&spacetime", "C_de_ii" ) );
+
+ rf.closeinput();
+
+ sprintf( fname, "%s/lpic.steps", p.file_path );
+
+ rf.openinput( fname );
+
+ cells_per_wl = atoi( rf.getinput( "spl" ) );
+ steps_per_period = atoi( rf.getinput( "spp" ) );
+
+ rf.closeinput();
+
+ bob.message("parameter read");
+
+ save(p);
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void input_spacetime::save( parameter &p )
+{
+ static error_handler bob("input_spacetime::save",errname);
+ ofstream outfile;
+
+ outfile.open(p.save_path_name,ios::app);
+
+ outfile << "spacetime" << endl;
+ outfile << "------------------------------------------------------------------" << endl;
+ outfile << "t_start : " << t_start << endl;
+ outfile << "t_stop : " << t_stop << endl;
+ outfile << "periods_t : " << periods_t << endl;
+ outfile << "steps_pp : " << steps_per_period << endl;
+ outfile << "x_start : " << x_start << endl;
+ outfile << "x_stop : " << x_stop << endl;
+ outfile << "periods_x : " << periods_x << endl;
+ outfile << "cells_per_wl : " << cells_per_wl << endl;
+ outfile << "x_offset : " << x_offset << endl;
+ outfile << "average : " << average << endl;
+ outfile << "size : " << size << endl;
+ outfile << "contour_1 : " << contour_1 << endl;
+ outfile << "contour_2 : " << contour_2 << endl;
+ outfile << "contour_3 : " << contour_3 << endl;
+
+ outfile.setf(ios::left);
+
+ outfile << "\n kw kt de di jx jy jz ex ey ez ed" <<endl;
+ outfile << "Q ";
+ outfile << ":" << setw(5) << Q_kw;
+ outfile << ":" << setw(5) << Q_kt;
+ outfile << ":" << setw(5) << Q_de;
+ outfile << ":" << setw(5) << Q_di;
+ outfile << ":" << setw(5) << Q_jx;
+ outfile << ":" << setw(5) << Q_jy;
+ outfile << ":" << setw(5) << Q_jz;
+ outfile << ":" << setw(5) << Q_ex;
+ outfile << ":" << setw(5) << Q_ey;
+ outfile << ":" << setw(5) << Q_ez;
+ // outfile << ":" << setw(5) << Q_bx;
+ outfile << ":" << setw(5) << Q_by;
+ outfile << ":" << setw(5) << Q_bz;
+ outfile << ":" << setw(5) << Q_edens;
+ outfile << ":" << setw(5) << Q_de_fi;
+ outfile << ":" << setw(5) << Q_de_ii << endl;
+
+ outfile << "C ";
+ outfile << ":" << setw(5) << C_kw;
+ outfile << ":" << setw(5) << C_kt;
+ outfile << ":" << setw(5) << C_de;
+ outfile << ":" << setw(5) << C_di;
+ outfile << ":" << setw(5) << C_jx;
+ outfile << ":" << setw(5) << C_jy;
+ outfile << ":" << setw(5) << C_jz;
+ outfile << ":" << setw(5) << C_ex;
+ outfile << ":" << setw(5) << C_ey;
+ outfile << ":" << setw(5) << C_ez;
+ // outfile << ":" << setw(5) << C_bx;
+ outfile << ":" << setw(5) << C_by;
+ outfile << ":" << setw(5) << C_bz;
+ outfile << ":" << setw(5) << C_edens;
+ outfile << ":" << setw(5) << C_de_fi;
+ outfile << ":" << setw(5) << C_de_ii << endl << endl;
+
+ outfile << "K_cut : " << K_cut << endl;
+ outfile << "W_cut : " << W_cut << endl << endl;
+
+ outfile.close();
+
+ bob.message("parameter written");
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::select( void )
+{
+ static error_handler bob("spacetime::select_plot",errname);
+
+ if (input.Q_de) {
+ xt_kt_kw( "de", input.C_de, 0, 0 );
+ write_idl_header_xt(input.C_de,0, "de", "spacetime-de", "!3n!De!N/n!Dc!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "de", "spacetime-kt-de", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "de", "spacetime-kw-de", "");
+ }
+
+ if (input.Q_di) {
+ xt_kt_kw( "di", input.C_di, 0, 0 );
+ write_idl_header_xt(input.C_di,0, "di", "spacetime-di", "!3n!Di!N/n!Dc!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "di", "spacetime-kt-di", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "di", "spacetime-kw-di", "");
+ }
+
+ if (input.Q_jx) {
+ xt_kt_kw( "jx", input.C_jx, 1, 0 );
+ write_idl_header_xt(input.C_jx,1, "jx", "spacetime-jx", "!3j!Dx!N/(e n!Dc!3 c)");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "jx", "spacetime-kt-jx", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "jx", "spacetime-kw-jx", "");
+ }
+
+ if (input.Q_jy) {
+ xt_kt_kw( "jy", input.C_jy, 1, 0 );
+ write_idl_header_xt(input.C_jy,1, "jy", "spacetime-jy", "!3j!Dy!N/(e n!Dc!3 c)");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "jy", "spacetime-kt-jy", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "jy", "spacetime-kw-jy", "");
+ }
+
+ if (input.Q_jz) {
+ xt_kt_kw( "jz", input.C_jz, 1, 0 );
+ write_idl_header_xt(input.C_jz,1, "jz", "spacetime-jz", "!3j!Dz!N/(e n!Dc!3 c)");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "jz", "spacetime-kt-jz", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "jz", "spacetime-kw-jz", "");
+ }
+
+ if (input.Q_ex) {
+ xt_kt_kw( "ex", input.C_ex, 1, 0 );
+ write_idl_header_xt(input.C_ex,1, "ex", "spacetime-ex", "!3E!Dx!N/E!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt,1, "ex", "spacetime-kt-ex", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0 , "ex", "spacetime-kw-ex", "");
+ }
+
+ if (input.Q_ey) {
+ xt_kt_kw( "ey", input.C_ey, 1, 0 );
+ write_idl_header_xt(input.C_ey,1, "ey", "spacetime-ey", "!3E!Dy!N/E!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "ey", "spacetime-kt-ey", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "ey", "spacetime-kw-ey", "");
+ }
+
+ if (input.Q_ez) {
+ xt_kt_kw( "ez", input.C_ez, 1, 0 );
+ write_idl_header_xt(input.C_ez,1, "ez", "spacetime-ez", "!3E!Dz!N/E!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "ez", "spacetime-kt-ez", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "ez", "spacetime-kw-ez", "");
+ }
+
+ if (input.Q_bx) {
+ xt_kt_kw( "bx", input.C_bx, 1, 0 );
+ write_idl_header_xt(input.C_bx,1, "bx", "spacetime-bx", "!3B!Dx!N/B!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "bx", "spacetime-kt-bx", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "bx", "spacetime-kw-bx", "");
+ }
+
+ if (input.Q_by) {
+ xt_kt_kw( "by", input.C_by, 1, 0 );
+ write_idl_header_xt(input.C_by,1, "by", "spacetime-by", "!3B!Dy!N/B!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "by", "spacetime-kt-by", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "by", "spacetime-kw-by", "");
+ }
+
+ if (input.Q_bz) {
+ xt_kt_kw( "bz", input.C_bz, 1, 0 );
+ write_idl_header_xt(input.C_bz,1, "bz", "spacetime-bz", "!3B!Dz!N/B!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "bz", "spacetime-kt-bz", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "bz", "spacetime-kw-bz", "");
+ }
+
+ if (input.Q_edens) {
+ xt_kt_kw( "edens", input.C_edens, 0, 0 );
+ write_idl_header_xt(input.C_edens,0, "edens", "spacetime-edens", "!3W!N/W!Dr!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "edens", "spacetime-kt-edens", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "edens", "spacetime-kw-edens", "");
+ }
+
+ if (input.Q_de_fi) {
+ xt_kt_kw( "de_fi", input.C_de_fi, 0, 1 );
+ write_idl_header_xt(input.C_de_fi,0, "de_fi", "spacetime-de_fi", "!3n!De!N/n!Dc!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "de_fi", "spacetime-kt-de_fi", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "de_fi", "spacetime-kw-de_fi", "");
+ }
+
+ if (input.Q_de_ii) {
+ xt_kt_kw( "de_ii", input.C_de_ii, 0, 1 );
+ write_idl_header_xt(input.C_de_ii,0, "de_ii", "spacetime-de_ii", "!3n!De!N/n!Dc!3");
+ if (input.Q_kt) write_idl_header_kt(input.C_kt, 1, "de_ii", "spacetime-kt-de_ii", "");
+ if (input.Q_kw) write_idl_header_kw(input.C_kw, 0, "de_ii", "spacetime-kw-de_ii", "");
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::xt_kt_kw( char *unit, float cut, int sign, int scale_write )
+{
+ static error_handler bob("spacetime::xt_kt_kw",errname);
+ char fname_xt[filename_size];
+ char fname_kt[filename_size];
+ char fname_kw[filename_size];
+
+ sprintf( fname_xt, "spacetime-%s", unit );
+ sprintf( fname_kt, "spacetime-kt-%s", unit );
+ sprintf( fname_kw, "spacetime-kw-%s", unit );
+
+ read_input_array_size( fname_xt ); // initializes x_steps_in, t_steps_in
+ // spp, spl, fnumber
+
+ printf( "reading %s ...",fname_xt ); fflush(stdout);
+ matrix_read = fmatrix( 0, t_steps_in-1, 0, x_steps_in-1 ); // array dimensions
+ read( fname_xt );
+ printf( "done\n" );
+
+ printf( "smoothing ... " ); fflush(stdout);
+ smooth( matrix_read );
+ printf( "done\n" );
+
+ printf( "scaling to byte ... " ); fflush(stdout);
+ matrix_write = ucmatrix( 0, t_steps_in-1, 0, x_steps_in-1 );
+ scale( cut, sign, matrix_read, matrix_write, t_steps_in, x_steps_in );
+ printf( "done\n" );
+
+ printf( "writing to disk ... " ); fflush(stdout);
+ write( fname_xt, scale_write );
+ bob.message("check point ktkw");
+ delete_ucmatrix( matrix_write, 0, t_steps_in-1, 0, x_steps_in-1 );
+ printf( "done\n" );
+
+ if (input.Q_kt) {
+
+ printf( "transforming x,t -> k,t ... " ); fflush(stdout);
+ kspace = fmatrix( 0, t_steps_in-1, 0, ft.steps_half-1 );
+ transform_k( matrix_read );
+ printf( "done\n" );
+
+ printf( "scaling to byte ... " ); fflush(stdout);
+ matrix_write = ucmatrix( 0, t_steps_in-1, 0, ft.steps_half-1 );
+ scale( input.C_kt, 0, kspace, matrix_write, t_steps_in, ft.steps_half );
+ delete_fmatrix( kspace, 0, t_steps_in-1, 0, ft.steps_half-1 );
+ printf( "done\n" );
+
+ printf( "writing to disk ... " ); fflush(stdout);
+ write_transform_k( fname_kt, matrix_write );
+ delete_ucmatrix( matrix_write, 0, t_steps_in-1, 0, x_steps_in-1 );
+ printf( "done\n" );
+ }
+ if (input.Q_kw) {
+
+ printf( "transforming x,t -> k,w ... " ); fflush(stdout);
+ transform_kw( matrix_read );
+ printf( "done\n" );
+
+ printf( "scaling to byte ... " ); fflush(stdout);
+ matrix_write = ucmatrix( 0, ft2d.steps_half_1-1, 0, ft2d.steps_2-1 );
+ scale( input.C_kw, 0, ft2d.power, matrix_write, ft2d.steps_half_1, ft2d.steps_2 );
+ printf( "done\n" );
+
+ printf( "writing to disk ... " ); fflush(stdout);
+ write_transform_kw( fname_kw, matrix_write );
+ delete_ucmatrix( matrix_write, 0, ft2d.steps_half_1-1, 0, ft2d.steps_2-1);
+ printf( "done\n" );
+ }
+ delete_fmatrix( matrix_read, 0, t_steps_in-1, 0, x_steps_in-1 );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::read_input_array_size( char *unit )
+{
+ static error_handler bob("spacetime::read_input_array_size",errname);
+
+ FILE *file;
+ char fname[ filename_size ];
+ float x_start, x_stop;
+ int x_steps;
+ int period;
+
+ fnumber = 0;
+ x_steps_in = 0;
+ do // read all spacetime file headers in order to determine
+ { // the dimension of the input array
+ sprintf( fname, "%s/%s-%d-%d", input_path, unit, fnumber+1, input.t_start );
+ file = fopen( fname, "rb" );
+ bob.message( "filename = ",fname);
+ if (file) {
+ fnumber++;
+
+ fread( &period, sizeof(int), 1, file );
+ fread( &spp, sizeof(int), 1, file );
+ fread( &x_start, sizeof(float), 1, file );
+ fread( &x_stop, sizeof(float), 1, file );
+ fread( &x_steps, sizeof(int), 1, file );
+
+ fclose( file );
+
+ if (x_steps_in==0 && x_steps>0) x_start_in = 1e-6 * floor( 1e6 * x_start );
+
+ if (x_steps>0 ) {
+ spl = (int) floor( (float) x_steps / (x_stop - x_start) + 0.5 );
+ x_stop_in = 1e-6 * floor( 1e6 * x_stop );
+ }
+
+ x_steps_in += x_steps;
+ }
+ }
+ while( file );
+
+ if (fnumber==0) bob.error( "no spacetime files found at time", input.t_start );
+
+ bob.message( "found", fnumber, "domain-spacetime file(s) at time", input.t_start );
+ bob.message( "spp =", spp );
+ bob.message( "x_start =", x_start_in );
+ bob.message( "x_stop =", x_stop_in );
+ bob.message( "x_steps =", x_steps_in );
+ bob.message( "spl =", spl );
+
+ if ( input.x_start < x_start_in ) input.x_start = x_start_in;
+ if ( input.x_stop > x_stop_in ) input.x_stop = x_stop_in;
+
+ t_steps_in = (input.t_stop - input.t_start) * spp;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::read( char *unit )
+{
+ static error_handler bob("spacetime::read",errname);
+
+ FILE *file;
+ char fname[ filename_size ];
+ float x_start, x_stop;
+ int x_steps;
+ int ti, fi, *x_steps_previous;
+ int period, period_in, spp_in;
+
+ x_steps_previous = new int [ t_steps_in ];
+
+ for( ti=0; ti<t_steps_in; ti++ ) x_steps_previous[ti] = 0;
+
+ for( period=input.t_start; period<input.t_stop; period++ )
+ {
+ bob.message("now at time ", period);
+ for( fi=1; fi<=fnumber; fi++ ) // now read the files
+ {
+ sprintf( fname, "%s/%s-%d-%d", input_path, unit, fi, period );
+ file = fopen( fname, "rb" );
+ if (!file) bob.error( "Cannot open file", fname );
+ else bob.message( "reading file", fname );
+
+ fread( &period_in, sizeof(int), 1, file );
+ fread( &spp_in, sizeof(int), 1, file );
+
+ if (period!=period_in) bob.error( "wrong period in file", fname );
+ if (spp!=spp_in) bob.error( "wrong number of time steps in file", fname );
+
+ for( ti=(period-input.t_start)*spp; ti<(period-input.t_start)*spp+spp; ti++ ) {
+ fread( &x_start, sizeof(float), 1, file );
+ fread( &x_stop, sizeof(float), 1, file );
+ fread( &x_steps, sizeof(int), 1, file );
+ fread( matrix_read[ti] + x_steps_previous[ti], sizeof(float), x_steps, file );
+ (x_steps_previous[ti]) += x_steps;
+ }
+
+ fclose( file );
+ }
+ }
+
+ bob.message( "all spacetime data read!" );
+
+ delete x_steps_previous;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::scale( float cut, int sign,
+ float **m, unsigned char **mw, int nt, int nx )
+ // sign = 1 : signed data
+ // sign = 0 : unsigned data ( de, di, edens )
+{
+ static error_handler bob("spacetime::scale",errname);
+
+ int INTMAX = 255;
+ int INTMAX_HALF = 127;
+ int overflow=0;
+ float data;
+ int ti, xi;
+
+ bob.message( "cut =", cut );
+
+ for( ti=0; ti<nt; ti++ ) { // scale and
+ for( xi=0; xi<nx; xi++ ) { // convert to unsigned char
+ if (sign==0) data = (float) fabs( m[ti][xi]/cut * INTMAX );
+ else data = (float) (INTMAX_HALF + INTMAX_HALF * m[ti][xi]/cut);
+ if (data>INTMAX) { data=INTMAX; overflow++; }
+ if (data<0) { data=0; }
+ mw[ti][xi] = (unsigned char) floor( data );
+ }
+ }
+
+ bob.message( "overflow =", (float)overflow/(nt*nx) );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::smooth( float **field )
+{
+ static error_handler bob("spacetime::smooth",errname);
+
+ float *s;
+ int offset = (int) floor( 0.5*( input.average - 1 ) + 0.1 );
+ int ti, xi, xj;
+
+ if ( input.average<3 ) return;
+
+ s = new float [x_steps_in];
+
+ for( ti=0; ti<t_steps_in; ti++ ) { // forall times
+
+ for( xi=0; xi<x_steps_in; xi++ ) s[xi] = field[ti][xi];
+
+ for( xi=offset; xi<=x_steps_in-offset; xi++ ) {
+
+ field[ti][xi] = 0;
+ for( xj=xi-offset; xj<=xi+offset; xj++ ) field[ti][xi] += s[xj] / input.average;
+ }
+ }
+
+ delete s;
+
+ bob.message( "smooth =", input.average );
+ bob.message( "smoothing done" );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write( char *unit, int scale_write )
+{
+ static error_handler bob("spacetime::write",errname);
+
+ FILE *file;
+ char filename[50];
+ float t;
+ float x, xmin, xmax;
+ int tu, to, xu, xo;
+ float data;
+ int i, j;
+ int k, l, flag, max_flag=0;
+
+ if ( input.x_start < x_start_in ) input.x_start = x_start_in;
+ if ( input.x_stop > x_stop_in ) input.x_stop = x_stop_in;
+
+ xmin = floor( (input.x_start - x_start_in) * spl + 0.5 );
+ xmax = floor( (input.x_stop - x_start_in) * spl + 0.5 );
+
+ bob.message( "t_start_out =", input.t_start );
+ bob.message( "t_stop_out =", input.t_stop );
+ bob.message( "x_start_out =", input.x_start );
+ bob.message( "x_stop_out =", input.x_stop );
+ bob.message( "spp =", spp );
+ bob.message( "spl =", spl );
+ bob.message( "t_steps_in =", t_steps_in );
+ bob.message( "xmin =", xmin );
+ bob.message( "xmax =", xmax );
+
+ sprintf( filename, "%s/%s", output_path, unit );
+ file = fopen( filename, "wb" );
+ if (!file) bob.error( "cannot open file", filename );
+
+ bob.message("check point 1");
+ vector_write = new unsigned char [input.size];
+
+ if (scale_write == 1) // de_fi and de_ii
+ {
+ bob.message("check point 2");
+
+ for( i=0; i<input.size; i++ ) // determine max_flag
+ {
+ for( j=0; j<input.size; j++ )
+ {
+ t = t_steps_in * i/input.size; if (t<1) t=1.001;
+ x = xmin + (xmax-xmin) * j/input.size; if (x<1) x=1.001;
+ to = (int) ceil( t ); tu = (int) floor( t );
+ xo = (int) ceil( x ); xu = (int) floor( x );
+ t = (float) to - t;
+ x = (float) xo - x;
+
+ flag = 0;
+ for( k = tu - int(ceil(t_steps_in/input.size));
+ k <= to + int(ceil(t_steps_in/input.size)); k++)
+ {
+ for( l = xu - int(ceil((xmax-xmin)/input.size));
+ l <= xo + int(ceil((xmax-xmin)/input.size)); l++)
+ {
+ if ( k >= 0 && k < t_steps_in && l >= 0 && l < x_steps_in &&
+ matrix_write[k][l]>0) flag ++;
+ }
+ }
+ max_flag = flag > max_flag ? flag:max_flag;
+ }
+ }
+ for( i=0; i<input.size; i++ ) // interpolate to output array dimensions
+ {
+ for( j=0; j<input.size; j++ )
+ {
+ t = t_steps_in * i/input.size; if (t<1) t=1.001;
+ x = xmin + (xmax-xmin) * j/input.size; if (x<1) x=1.001;
+ to = (int) ceil( t ); tu = (int) floor( t );
+ xo = (int) ceil( x ); xu = (int) floor( x );
+ t = (float) to - t;
+ x = (float) xo - x;
+
+ flag = 0;
+ for( k = tu - int(ceil(t_steps_in/input.size));
+ k <= to + int(ceil(t_steps_in/input.size)); k++)
+ {
+ for( l = xu - int(ceil((xmax-xmin)/input.size));
+ l <= xo + int(ceil((xmax-xmin)/input.size)); l++)
+ {
+ if ( k >= 0 && k < t_steps_in && l >= 0 && l < x_steps_in &&
+ matrix_write[k][l]>0) flag ++;
+ }
+ }
+ if (flag>0) vector_write[j] = (unsigned char) floor(255*flag/max_flag+0.5);
+ else vector_write[j] = (unsigned char) 0;
+ }
+ fwrite( vector_write, sizeof(unsigned char), input.size, file );
+ }
+ }
+
+ if (scale_write == 0) // all the others
+ {
+ bob.message("check point 3");
+
+ for( i=0; i<input.size; i++ ) // interpolate to output array dimensions
+ {
+ for( j=0; j<input.size; j++ )
+ {
+ t = t_steps_in * i/input.size; if (t<1) t=1.001;
+ x = xmin + (xmax-xmin) * j/input.size; if (x<1) x=1.001;
+ to = (int) ceil( t ); tu = (int) floor( t );
+ xo = (int) ceil( x ); xu = (int) floor( x );
+ t = (float) to - t;
+ x = (float) xo - x;
+ data = t*x*matrix_write[tu][xu] +
+ t*(1.0-x)*matrix_write[tu][xo];
+ data += (1.0-t)*x*matrix_write[to][xu] +
+ (1.0-t)*(1.0-x)*matrix_write[to][xo];
+ vector_write[j] = (unsigned char) floor( data + 0.5 );
+ }
+ fwrite( vector_write, sizeof(unsigned char), input.size, file );
+ }
+ }
+
+ fclose( file );
+
+ bob.message("check point 4");
+ delete vector_write;
+ bob.message("check point 5");
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write_idl_header_xt( float cut, int sign, char *idlname, char *unit,
+ char *axislabel )
+{
+ static error_handler bob("spacetime::write_idl_header",errname);
+
+ FILE *file;
+ int i, j;
+ unsigned char low;
+
+ char fname[filename_size];
+ sprintf( fname, "%s/scale-idl", output_path );
+ file = fopen( fname, "wb" );
+
+ for( i=0; i<=255; i++ ) // write color table
+ {
+ low = (unsigned char) i;
+ for( j=1; j<=30; j++ ) fwrite( &(low), sizeof(unsigned char), 1, file );
+ }
+ fclose( file );
+
+ sprintf( fname, "%s/idl_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl_%s", idlname );
+ fprintf( file, "\n\n xoffset = %.4f", input.x_offset );
+ fprintf( file, "\n\n xmin = %.2f - xoffset", input.x_start );
+ fprintf( file, "\n xmax = %.2f - xoffset", input.x_stop );
+ fprintf( file, "\n xname = \"!3x/!7k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %d", input.t_start );
+ fprintf( file, "\n ymax = %d", input.t_stop );
+ fprintf( file, "\n yname = \"!3t/!7s!3\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 4" );
+ fprintf( file, "\n\n unit = \"%s\"\n\n", unit );
+ fclose( file );
+
+ sprintf( fname, "%s/idl2ps_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl2ps_%s", idlname );
+ fprintf( file, "\n\n xoffset = %.4f", input.x_offset );
+ fprintf( file, "\n\n xmin = %.2f - xoffset", input.x_start );
+ fprintf( file, "\n xmax = %.2f - xoffset", input.x_stop );
+ fprintf( file, "\n xname = \"!3x/!7k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %d", input.t_start );
+ fprintf( file, "\n ymax = %d", input.t_stop );
+ fprintf( file, "\n yname = \"!3t/!7s!3\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 0" );
+ fprintf( file, "\n\n unit = \"%s\"", unit );
+ fprintf( file, "\n\n idl2eps = \"%s.eps\"\n\n", unit );
+ fclose( file );
+
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::transform_k( float **matrix_read )
+{
+ static error_handler bob("spacetime::transform",errname);
+
+ int i, j;
+ int steps = input.periods_x * input.cells_per_wl;
+
+ bob.message( "steps =", (double)steps, " steps_ft =", (double)ft.steps );
+
+ for( i=0; i<t_steps_in; i++ ) {
+ ft.RealFt( matrix_read[i] );
+ for( j=0; j<ft.steps_half; j++ ) {
+ kspace[i][j] = ft.power[j];
+ // kspace[i][j] = ft.co[j];
+ // kspace[i][j] = ft.si[j];
+ }
+ }
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write_transform_k( char *unit, unsigned char **m )
+{
+ static error_handler bob("spacetime::write_transform_k",errname);
+
+ FILE *file;
+ char filename[50];
+ float t;
+ float x, xmin, xmax, kmax;
+ int tu, to, xu, xo;
+ float data;
+ int i, j;
+ int nx=ft.steps_half;
+
+ kmax = ft.df * nx;
+ input.K_cut = (int) fabs( input.K_cut );
+
+ xmin = 0;
+
+ if (input.K_cut==0) { xmax = nx; input.K_cut = kmax; }
+ else if (input.K_cut<=kmax) { xmax = (int) floor( 1.0*nx*input.K_cut/kmax ); }
+ else { xmax = nx; input.K_cut = kmax; }
+
+ sprintf( filename, "%s/%s", output_path, unit );
+ file = fopen( filename, "wb" );
+ if (!file) bob.error( "cannot open file", filename );
+
+ vector_write = new unsigned char [input.size];
+
+ for( i=0; i<input.size; i++ ) // interpolate to output array dimensions
+ {
+ for( j=0; j<input.size; j++ )
+ {
+ t = t_steps_in * i/input.size; if (t<1) t=1.001;
+ x = xmin + (xmax-xmin) * j/input.size; if (x<1) x=1.001;
+ to = (int) ceil( t ); tu = (int) floor( t );
+ xo = (int) ceil( x ); xu = (int) floor( x );
+ t = (float) to - t;
+ x = (float) xo - x;
+ data = t*x*m[tu][xu] + t*(1.0-x)*m[tu][xo];
+ data += (1.0-t)*x*m[to][xu] + (1.0-t)*(1.0-x)*m[to][xo];
+ vector_write[j] = (unsigned char) floor( data + 0.5 );
+ }
+ fwrite( vector_write, sizeof(unsigned char), input.size, file );
+ }
+
+ fclose( file );
+
+ delete vector_write;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write_idl_header_kt( float cut, int sign, char *idlname, char *unit,
+ char *axislabel )
+{
+ static error_handler bob("spacetime::write_idl_header_kt",errname);
+
+ FILE *file;
+ int i, j;
+ unsigned char low;
+
+ char fname[filename_size];
+ sprintf( fname, "%s/scale-idl", output_path );
+ file = fopen( fname, "wb" );
+
+ for( i=0; i<=255; i++ ) // write color table
+ {
+ low = (unsigned char) i;
+ for( j=1; j<=30; j++ ) fwrite( &(low), sizeof(unsigned char), 1, file );
+ }
+ fclose( file );
+
+ sprintf( fname, "%s/idl_kt_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl_kt_%s", idlname );
+ fprintf( file, "\n\n xmin = %.2f", 0.0 );
+ fprintf( file, "\n xmax = %.2f", input.K_cut );
+ fprintf( file, "\n xname = \"!3k/!3k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %d", input.t_start );
+ fprintf( file, "\n ymax = %d", input.t_stop );
+ fprintf( file, "\n yname = \"!3t/!7s!3\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 4" );
+ fprintf( file, "\n\n unit = \"%s\"\n\n", unit );
+ fclose( file );
+
+ sprintf( fname, "%s/idl2ps_kt_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl2ps_kt_%s", idlname );
+ fprintf( file, "\n\n xmin = %.2f", 0.0 );
+ fprintf( file, "\n xmax = %.2f", input.K_cut );
+ fprintf( file, "\n xname = \"!3k/!3k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %d", input.t_start );
+ fprintf( file, "\n ymax = %d", input.t_stop );
+ fprintf( file, "\n yname = \"!3t/!7s!3\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 0" );
+ fprintf( file, "\n\n unit = \"%s\"", unit );
+ fprintf( file, "\n\n idl2eps = \"%s.eps\"\n\n", unit );
+ fclose( file );
+
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::transform_kw( float **matrix_read )
+{
+ static error_handler bob("spacetime::transform_kw",errname);
+
+ int i, j;
+ int steps_x = input.periods_x * input.cells_per_wl;
+ int steps_t = input.periods_t * input.steps_per_period;
+
+ bob.message( "steps =", (double)steps_x, " steps_ft =", (double)ft2d.steps_2 );
+ bob.message( "steps_t", (double)steps_t, " steps_ft =", (double)ft2d.steps_1 );
+
+ ft2d.RealFt( matrix_read );
+
+ FILE *f;
+ f=fopen( "ft2d.dat", "w" );
+ for( i=0; i<ft2d.steps_half_1; i+= (int) floor(0.1 * ft2d.steps_half_1) ) {
+ fprintf( f, "\n" );
+ for( j=0; j<ft2d.steps_half_2; j+=10 ) {
+ fprintf( f, "\n %.4e %.3e", ft2d.frequency_2[j], ft2d.power[i][j] );
+ }
+ }
+ fclose( f );
+ // result in ft2d.power[i][j], ft2d.co[i][j], ft2d.si[i][j]
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write_transform_kw( char *unit, unsigned char **m )
+{
+ static error_handler bob("spacetime::write_transform_kw",errname);
+
+ FILE *file;
+ char filename[50];
+ float t, tmin, tmax, wmax;
+ float x, xmin, xmax, kmax;
+ int tu, to, xu, xo;
+ float data;
+ int i, j;
+
+ wmax = ft2d.df_1 * ft2d.steps_half_1;
+ kmax = ft2d.df_2 * ft2d.steps_half_2;
+
+ printf( "wmax = %.3f\n", wmax );
+ printf( "kmax = %.3f\n", kmax );
+
+ if (input.K_cut>0 && input.K_cut<=kmax) {
+ xmin = ft2d.steps_half_2 - (int) floor(1.0*input.K_cut/ft2d.df_2);
+ xmax = ft2d.steps_half_2 + (int) floor(1.0*input.K_cut/ft2d.df_2);
+ }
+ else {
+ xmin = 0;
+ xmax = ft2d.steps_2;
+ }
+
+ if (input.W_cut>0 && input.W_cut<=wmax) {
+ tmin = 0;
+ tmax = (int) floor( 1.0*input.W_cut/ft2d.df_1 );
+ }
+ else {
+ tmin = 0;
+ tmax = ft2d.steps_half_1;
+ }
+
+ sprintf( filename, "%s/%s", output_path, unit );
+ file = fopen( filename, "wb" );
+ if (!file) bob.error( "cannot open file", filename );
+
+ vector_write = new unsigned char [input.size];
+
+ for( i=0; i<input.size; i++ ) // interpolate to output array dimensions
+ {
+ for( j=0; j<input.size; j++ )
+ {
+ t = tmin + (tmax-tmin) * i/input.size; if (t<1) t=1.001;
+ x = xmin + (xmax-xmin) * j/input.size; if (x<1) x=1.001;
+ to = (int) ceil( t ); tu = (int) floor( t );
+ xo = (int) ceil( x ); xu = (int) floor( x );
+ t = (float) to - t;
+ x = (float) xo - x;
+ data = t*x*m[tu][xu] + t*(1.0-x)*m[tu][xo];
+ data += (1.0-t)*x*m[to][xu] + (1.0-t)*(1.0-x)*m[to][xo];
+ vector_write[j] = (unsigned char) floor( data + 0.5 );
+ }
+ fwrite( vector_write, sizeof(unsigned char), input.size, file );
+ }
+
+ fclose( file );
+
+ delete vector_write;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+
+void spacetime::write_idl_header_kw( float cut, int sign, char *idlname, char *unit,
+ char *axislabel )
+{
+ static error_handler bob("spacetime::write_idl_header_kw",errname);
+
+ FILE *file;
+ int i, j;
+ unsigned char low;
+
+ char fname[filename_size];
+ sprintf( fname, "%s/scale-idl", output_path );
+ file = fopen( fname, "wb" );
+
+ for( i=0; i<=255; i++ ) // write color table
+ {
+ low = (unsigned char) i;
+ for( j=1; j<=30; j++ ) fwrite( &(low), sizeof(unsigned char), 1, file );
+ }
+ fclose( file );
+
+ sprintf( fname, "%s/idl_kw_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl_kw_%s", idlname );
+ fprintf( file, "\n\n xmin = %.2f", -input.K_cut );
+ fprintf( file, "\n xmax = %.2f", input.K_cut );
+ fprintf( file, "\n xname = \"!3k/!3k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %.3f", 0.0 );
+ fprintf( file, "\n ymax = %.3f", input.W_cut );
+ fprintf( file, "\n yname = \"!7x/!7x!3!D0\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 4" );
+ fprintf( file, "\n\n unit = \"%s\"\n\n", unit );
+ fclose( file );
+
+ sprintf( fname, "%s/idl2ps_kw_%s.header", output_path, idlname );
+ file = fopen( fname, "w" );
+ fprintf( file, "pro idl2ps_kw_%s", idlname );
+ fprintf( file, "\n\n xmin = %.2f", -input.K_cut );
+ fprintf( file, "\n xmax = %.2f", input.K_cut );
+ fprintf( file, "\n xname = \"!3k/!3k!3!D0!3\"" );
+ fprintf( file, "\n ymin = %.3f", 0.0 );
+ fprintf( file, "\n ymax = %.3f", input.W_cut );
+ fprintf( file, "\n yname = \"!7x/!7x!3!D0\"" );
+ fprintf( file, "\n zmax = %.2e", cut );
+ if (sign==0) fprintf( file, "\n zmin = %.2e", 0.0 );
+ else fprintf( file, "\n zmin = %.2e", -cut );
+ fprintf( file, "\n zname = \"%s\"", axislabel );
+ fprintf( file, "\n\n level1 = %.2e", input.contour_1 );
+ fprintf( file, "\n level2 = %.2e", input.contour_2 );
+ fprintf( file, "\n level3 = %.2e", input.contour_3 );
+ fprintf( file, "\n\n dim = %d", input.size );
+ fprintf( file, "\n\n color = 0" );
+ fprintf( file, "\n\n unit = \"%s\"", unit );
+ fprintf( file, "\n\n idl2eps = \"%s.eps\"\n\n", unit );
+ fclose( file );
+
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+//eof
+
+
+
+
+
+
+
+
+
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-05-22 08:35:00 +0000