SOS: gnuplot fdtd的一个问题求助 perl vs python
我用perl和python写了相同功能的一段程序，计算一维fdtd,用gnuplot动态显示，可是python的数据没有显示出来，看横纵坐标的变化数据是正确收到了的，如最后的图片，求大神指点，谢谢。
1 #!/usr/bin/perl 2 3 # Author : Leon Email: yangli0534@gmail.com 4 # fdtd simulation , plotting with gnuplot, writting in perl 5 # perl and gnuplot software packages should be installed before running this program 6 # 1d fdtd with absorbing boundary and TFSF boundary between [49] and [50] 7 # lossy dielectric material localted at > ez[150] 8 #use Time::HiRes qw(sleep); 9 #use autodie qw(:all); 10 use strict; 11 use warnings; 12 13 #print "\@\n"; 14 #my $terminal = ""; 15 # open GNUPLOT_TERM, "echo 'show terminal;' | gnuplot 2>&1 |"; 16 # while (<GNUPLOT_TERM>) { 17 # if (m/terminal type is (\w+)/) { 18 # $terminal=$1; 19 # } 20 # } 21 # close GNUPLOT_TERM; 22 23 # unfortunately, the wxt terminal type does not support positioning. 24 # hardcode it... 25 # $terminal = "x11"; 26 27 open my $PIPE ,"| gnuplot " || die "Can't initialize gnuplot number \n"; 28 29 print $PIPE "set size 0.85, 0.85\n"; 30 print $PIPE "set term png size 600, 400\n"; 31 32 my $title = "fdtd simulation by leon,yangli0534\\\\@"."gmail.com"; 33 print $PIPE "set terminal gif animate\n";# terminal type: png 34 print $PIPE "set output \"fdtd_simulation_abc_tfsf_match_diel.gif\"\n";#output file name 35 print $PIPE "set title \"{/Times:Italic $title}\"\n";# title name and font 36 #print $PIPE "set title \"fdtd simulation by leon,yangli0534\\\\@ gmail.com\"\n";# title name and font 37 print $PIPE "set title font \",10\" norotate tc rgb \"white\"\n"; 38 print $PIPE "unset key\n"; 39 print $PIPE "set tics textcolor rgb \"white\"\n";# text color 40 print $PIPE "set border lc rgb \"orange\"\n"; 41 print $PIPE "set grid lc rgb\"orange\"\n"; 42 print $PIPE "set object 1 rectangle from screen 0,0 to screen 1,1 fc rgb \"gray10\" behind\n";#background color 43 print $PIPE "set xlabel\" {/Times:Italic distance: wave length}\" tc rgb \"white\" \n";# xlabel 44 print $PIPE "set ylabel\"{/Times:Italic amplitude: v}\" tc rgb \"white\"\n";#ylabel 45 print $PIPE "set autoscale\n"; 46 47 my $size = 400;#physical distance 48 my @ez;#electric field 49 my @hy;#magnetic field 50 my @ceze;# 51 my @cezh;# 52 my @chye;# 53 my @chyh;# 54 my $LOSS_LAYER = 250; 55 56 #my @epsR; 57 my $LOSS = 0.01; 58 my $imp0 = 377.0; 59 #initalization 60 for (my $i = 0; $i < $size; $i++){ 61 $ez[$i] = 0.0; 62 $hy[$i] = 0.0; 63 if ($i < 100) { 64 #$epsR[$i] = 1.0; 65 $ceze[$i] = 1.0; 66 $cezh[$i] = $imp0; 67 } 68 elsif($i < $LOSS_LAYER){ 69 #$epsR[$i] = 1.0; 70 $ceze[$i] = 1.0; 71 $cezh[$i] = $imp0/9.0; 72 } 73 else { 74 #$epsR[$i] = 9.0; 75 $ceze[$i] = (1.0-$LOSS)/(1.0+$LOSS); 76 $cezh[$i] = $imp0 / 9 /(1.0+$LOSS); 77 } 78 } 79 80 for ( my $i = 0; $i < $size; $i++){ 81 82 83 if($i < $LOSS_LAYER){ 84 #$epsR[$i] = 1.0; 85 $chye[$i] = 1.0/$imp0; 86 $chyh[$i] = 1.0; 87 } 88 else { 89 #$epsR[$i] = 9.0; 90 $chye[$i] = 1.0/ $imp0 /(1.0+$LOSS); 91 $chyh[$i] = (1.0-$LOSS)/(1.0+$LOSS); 92 } 93 } 94 my $qTime; 95 my $MaxTime = 18050; 96 my $pi = 3.141592563589793; 97 print $PIPE "set xrange [0:$size-1]\n"; 98 my $mm = 0; 99 100 #do time stepping 101 for($qTime = 0; $qTime < $MaxTime; $qTime+=1){ 102 103 # update magnetic field 104 #$hy[$size - 1] = $hy[$size - 2 ];#abc 105 for( $mm = 0; $mm < $size - 1; $mm++){ 106 $hy[$mm] = $hy[$mm]*$chyh[$mm] + ($ez[$mm+1] - $ez[$mm])*$chye[$mm]; 107 } 108 $hy[49] -=exp(-($qTime - 30.0)*($qTime - 30.0)/100.0)/$imp0; 109 # update electric field 110 $ez[0] = $ez[1];#abc 111 #$ez[$size-1] = $ez[$size-2]; 112 for( $mm = 1; $mm < $size-1 ; $mm++){ 113 $ez[$mm] = $ez[$mm]*$ceze[$mm] + ($hy[$mm] - $hy[$mm-1])*$cezh[$mm]; 114 } 115 116 if($qTime % 20 == 0){ 117 118 print $PIPE "plot \"-\" w l lw 1.5 lc rgb \"green\"\n"; 119 my $cnt = 0; 120 for my $elem ( @ez) { 121 #print " ".$elem; 122 print $PIPE $cnt." ".$elem."\n"; 123 $cnt += 1; 124 } 125 print $PIPE "e\n"; 126 } 127 #hardwire a source 128 $ez[50] += exp(-($qTime +0.5-(-0.5)- 30.0)*($qTime +0.5-(-0.5)- 30.0)/100.0); 129 } 130 131 #print $PIPE "set terminal x11\n"; 132 133 print $PIPE "set output\n"; 134 135 close($PIPE);
正常结果
python的程序
1 #!/usr/bin/env python 2 3 import sys 4 import math 5 import os 6 # Author : Leon Email: yangli0534@gmail.com 7 # fdtd simulation , plotting with gnuplot, writting in python 8 # perl and gnuplot software packages should be installed before running this program 9 # 1d fdtd with absorbing boundary and TFSF boundary between [49] and [50] 10 # lossy dielectric material localted at > ez[150] 11 gp = os.popen('gnuplot','w') 12 #gp .write('plot sin(x)\n'); 13 #gp .flush(); 14 #p .close(); 15 gp.write('set size 0.85, 0.85\n') 16 gp.write('set term png size 600, 400\n') 17 title = 'fdtd simulation by leon,yangli0534\\\\@gmail.com' 18 gp.write('set terminal gif animate\n') 19 gp.write('set output "fdtd_simulation_abc_tfsf_match_diel_py.gif"\n') 20 gp.write(''.join(['set title "{/Times:Italic ',title, '}"\n'])) 21 gp.write('set title font ",10" norotate tc rgb "white"\n') 22 gp.write('unset key\n') 23 gp.write('set tics textcolor rgb "white"\n') 24 gp.write('set border lc rgb "orange"\n') 25 gp.write('set grid lc rgb"orange"\n') 26 gp.write('set object 1 rectangle from screen 0,0 to screen 1,1 fc rgb "gray10" behind\n') 27 gp.write('set xlabel" {/Times:Italic distance: wave length}" tc rgb "white" \n') 28 gp.write('set ylabel"{/Times:Italic amplitude: v}" tc rgb "white"\n') 29 gp.write('set autoscale\n'); 30 31 size = 400#physical distance 32 ez=size * [0.00]#electric field 33 hy=size * [0.00]#magnetic field 34 ceze=size * [0.00]# 35 cezh=size * [0.00]# 36 chye=size * [0.00]# 37 chyh=size * [0.00]# 38 imp0 = 377.00 39 LOSS = 0.01 40 LOSS_LAYER = 250 41 MaxTime = 18000 42 cnt = 0 43 elem = 0.00000 44 gp.write(''.join(['set xrange [0:',str(size),'-1]\n'])); 45 for i in range(0,size): 46 ez[i] = 0.0 47 hy[i] = 0.0 48 if (i < 100): 49 #$epsR[$i] = 1.0; 50 ceze[i] = 1.0 51 cezh[i] = imp0 52 elif(i < LOSS_LAYER): 53 #$epsR[$i] = 1.0; 54 ceze[i] = 1.0 55 cezh[i] = imp0/9.0 56 else : 57 #$epsR[$i] = 9.0; 58 ceze[i] = (1.0-LOSS)/(1.0+LOSS) 59 cezh[i] = imp0 / 9 /(1.0+LOSS) 60 if( i < LOSS_LAYER): 61 chye[i] = 1.0/imp0 62 chyh[i] = 1.0 63 else: 64 chye[i] = 1.0/imp0/(1.0+LOSS) 65 chyh[i] = (1.0-LOSS)/(1.0+LOSS) 66 for qTime in range(0, MaxTime): 67 # update magnetic field 68 for mm in range(0, size-1): 69 hy[mm] = hy[mm]*chyh[mm] + (ez[mm+1]-ez[mm])*chye[mm] 70 hy[49] = hy[49]-math.exp(-(qTime - 30.0)*(qTime - 30.0)/100.0)/imp0 71 # update electric field 72 ez[0] = ez[1]#abc 73 #$ez[$size-1] = $ez[$size-2]; 74 for mm in range(1, size-1): 75 ez[mm] = ez[mm]*ceze[mm] + (hy[mm] - hy[mm-1])*cezh[mm] 76 if(qTime % 20 == 0): 77 gp.write('plot "-" w l lw 1.5 lc rgb "white"\n') 78 cnt = 0 79 for elem in ez: 80 #gp.write([cnt,' ',elem,'\n']) 81 gp.write(''.join([str(cnt),' ',str(elem),'\n'])) 82 #gp.write(' ') 83 #gp.write(str(elem)) 84 #gp.write('\n') 85 cnt += 1 86 gp.write('e\n') 87 ez[50] = ez[50]+math.exp(-(qTime +0.5-(-0.5)- 30.0)*(qTime +0.5-(-0.5)- 30.0)/100.0); 88 gp.write('set output\n') 89 gp.close()
结果就是
OPTIMISM, PASSION & HARDWORK
