#!/opt/local/bin/perl use Math::Trig; #================================================= # perl script to demonstrate 'ioa' files from OpenGGCM. # Also plots fields and makes movies. # # Requirements: Perl5 # gnuplot version >=5, # ffmpeg # ps2pdf (psutils package) # pdfcrop # ImageMagick # Jimmy Raeder, June 2018 j.raeder@unh.edu # # ioa file structure: # for each variable: # # 1. 1024 character header # The header is all printable ascii # and contains metadata in the form 'var=value' # These variables are required in the header: # variable=XXXX variable name # nphi=nnn points in the phi (azimuthal) direction # 0-360 degrees, where 0 is midnight # ntheta= points in the theta (colatitude) direction # 0-180 degrees, where 0 is the north pole # All in SM coordinates # # There can be more metadata in the header, for example UTTIME, unit, etc. # The first header of a file can easily be viewed with # 'head -b 1024 filename | more' # After the metadata the header is filled with ascii zeros. # # # These variables are available as a minimum: # variable=pot potential in Volt # variable=prec_e_fe_1 electron precipitation flux, diffuse precipitation in W/m**2 # variable=prec_e_fe_2 electron precipitation flux, discrete precipitation in W/m**2 # variable=prec_e_e0_1 corresponding mean energy in eV # variable=prec_e_e0_2 corresponding mean energy in eV # variable=ftv_rcm flux tube volume. <0: open tubes, >0: closed tubes # # There can be more variables. # # 2. following the header come nphi*ntheta*4 bytes of binary (float/real*4) data # for the variable values on the spherical grid. # # The file design allows for concatenating files. # OpenGGCM writes one file for one time. # The file naming convention is: RUN.ioa.NNNNNN # where: # RUN: unique identifier of the simulation runa. # ioa: file type (there are also iof, ioc, 3df, px_xx and other types from OpenGGCM) # NNNNNN: seconds since the start of the run. The cadence can be set # in the OpenGGCM 'runme' file, right now it is 30 seconds. # The first 2 hours are basically the OpenGGCM startup (depends on the # event and purpose), so use 007200, 007230, .... # Look for UTTIME=.... in the header, given is ISO 8601 format. # # #================================================= # remote('target=jraeder@jr.sr.unh.edu rdir=tmp.remote'); #... run me on another machine #..... make list of files/times @a=glob("data/*ioa.02*");foreach(@a){ if(/^(.+)\.(\d{6})$/){ $R=$1; push(@S,$2); }} #...... loop over variables w/plot parameters foreach $w (qw/ pot:0:sym:-110:110:0.001 ftv_rcm:0:sym:0:0.5:1 prec_e_fe_1:1:log:8.0e-3:11.1:1000 /){ ($what,$log,$cmap,$min,$max,$scale)=split(':',$w); #...... make frames, save the frames $dd="dir.frames"; ssy("mkdir -p $dd"); ssy("/bin/rm -vf $dd/*"); foreach $s (@S){ $fn="$R.$s"; #.... filename print STDERR "next: ($what,$log,$cmap,$min,$max,$scale) |$fn|\n"; # exit; # @v=ioa_get($fn,'show'); # exit; #.... just show the headers $o="$dd/tmp_${what}_$s.jpg"; if( ! ( -e $o ) ){ #.... output frame #..... read the record, returns array w/header followed by all field values @v=ioa_get($fn,$what); $h=pop(@v); # print STDERR "\n\n $h \n\n"; #exit; # header $l="toplab=${fn}~$what~$s scale=$scale"; #..... top label, tildes are stripped #..... make the plot plot_frame("$h type=io-polar-N maxcolat=50 min=$min max=$max cmap=$cmap log=$log $l out=$o",@v); # exit ; }} make_movie("$dd/tmp_${what}_*jpg","$what.mp4"); #... make movie for this variable ssy("/bin/rm -vf tmp_*jpg tmpm*"); } #------------------------------------------------- sub plot_frame{ local($v,@r)=@_; #------------------------------------------------- local($n,$nphi,$nphi1,$ntheta,$ntheta1,$k,$levels,$min,$max); local($type,$variable,$zz,$it,$ip,$the,$phi,$r,$rmax,$nphi,$nphi1,$ntheta,$ntheta1); local($rlab,$levels,$cl); $toplab=parsev($v,'toplab=top-label'); $toplab=~s/\_/\\_/g; $toplab=~s/\~/ /g; $type=parsev($v,'type=latlonbox'); $rmax=parsev($v,'maxcolat=50.1')+0.1; $nphi=parsev($v,'nphi=10'); $ntheta=parsev($v,'ntheta=10'); $variable=parsev($v,'variable=none'); $UTTIME=parsev($v,'UTTIME=1900-01-01T00:00:00.00Z'); $nphi1=$nphi-1; $ntheta1=$ntheta-1; $min=parsev($v,'min=-1'); $max=parsev($v,'max=1'); $scale=parsev($v,'scale=1'); $dz=($max-$min)/200; $levels="${min},${dz},$max"; if($type eq 'io-polar-N'){ #...... make line styles for contours # by default gnuplot gives them all different colors, annoying $st=''; foreach (1..20){ $st .= "set style line $_ lw 0.7 lc rgb 'black' \n"; } #.... make tmp file w/data, convert to MLT and lat from colat open(FF,">tmp.1"); $z1=1e33; $z2=-$z1; $pi=4.0*atan(1.0);$rad=$pi/180; $np4=int($nphi/4); #... rotate for noon on top foreach $it (0..$ntheta1){ foreach $ip (0..$nphi1){ $phi=$rad*360*($ip/$nphi1); $ip=($ip+$np4)%$nphi; $k=$ip+$nphi*$it; $z=$r[$k]*$scale; $z1=min($z1,$z); $z2=max($z,$z2); $r=$it; if($r<$rmax){ $x=$r*cos($phi); $y=$r*sin($phi); #.. the color plot is actually Cart print FF "$x $y $z\n"; } } print FF "\n"; } close(FF); #..... make some labels $xlab=''; $z2=sprintf("%10.3g",$z2); $z1=sprintf("%10.3g",$z1); $xlab .= "set label '$z1' at screen 0.87,0.10 right font 'Helvetica,15' \n"; $xlab .= "set label '$z2' at screen 0.87,0.90 right font 'Helvetica,15' \n"; $_=$UTTIME; if(/^(.+)T(.+)$/){ $xlab .= "set label '$1' at screen 0.09,0.90 center font 'Helvetica,15' \n"; } $_=$UTTIME; if(/^(.+)T(.+)$/){ $xlab .= "set label '$2' at screen 0.09,0.88 center font 'Helvetica,15' \n"; } $rlab=$rmax*1.05; $dz=($max-$min)/12; if($log>0){ $dz=sqrt(10); } # #...... color maps, log $cl=''; #..... not yet happy with the sym color map if($cmap eq 'sym'){$cl .= "set palette model HSV functions 0.7+0.3*gray, 2*(abs(gray-0.5)**1.1), 1-1.0*abs(gray-0.5)**2\n"; } if($cmap eq 'log'){$cl .= "set palette rgbformulae 33,13,10\n"; } #.... standard rainbow if($log==1){ $cl .= "set logscale cb 10\n"; } #...... make gnuplot input $_=" set terminal postscript eps size 5,5 enhanced color font 'Helvetica,20' linewidth 3 set lmargin at screen 0.05; set rmargin at screen 0.85; set bmargin at screen 0.1; set tmargin at screen 0.9 set pm3d map unset key set label '$toplab' at screen 0.5,0.95 center font 'Helvetica,17' tc rgb 'black' $xlab set multiplot # plot the heatmap set isosamples 500 #...... for contour lines set cntrparam levels increment ${min},${dz},$max set contour surface #...... to make all contourlines same color $st set style increment user unset border; unset xtics; unset ytics; set angles degree; $cl set cbrange [${min}:${max}] set colorbox user origin 0.9,0.1 size 0.03,0.8 splot \'tmp.1\' using 1:2:3 lc rgb 'black' # now plot the polar grid only set style line 99 lc rgb 'black' lw 0.5 set grid polar front ls 99, ls 99 set polar set rrange[0:${rmax}] set rtics (10,20,30,40,50) scale 0 rotate by 90 font 'Helvetica,12' textcolor rgb 'grey' offset -0.8 unset parametric set for [i=0:330:30] label at first $rlab*cos(i), first $rlab*sin(i) center sprintf('%02d',((i+90)/15)%24) plot NaN w l #... plot labes only, no field unset multiplot # show label "; } if($type eq 'latlonbox'){ #..... this on is a bit ugly #.... make tmp file w/data, convert to MLT and lat from colat open(FF,">tmp.1"); foreach $it (0..$ntheta1){ foreach $ip (0..$nphi1){ $mlt=($ip/$nphi1)*24+12; $the=90-180*($it/$ntheta1); $k=$ip+$nphi*($ntheta1-$it); $zz=$r[$k]; print FF "$mlt $the $zz\n"; } print FF "\n"; } close(FF); $_=" set terminal postscript eps size 6,5 enhanced color font 'Helvetica,20' linewidth 3 #set lmargin at screen 0.10; set rmargin at screen 0.80; set bmargin at screen 0.1; set tmargin at screen 0.6 # set multiplot # # title does not show ???? set title \"$UTTIME $variable\" # set contour base; set cntrparam levels incremental $levels # unset surface; set table tmp_cont2.dat # splot 'tmp.1'; unset table; reset; #....... makes contour lines #set xrange [0:${nphi1}]; set yrange [0:${ntheta1}] set xrange [12:36]; set yrange [-90:90] unset key; set cbrange [${min}:${max}] set view map set palette rgbformulae 33,13,10 set tics scale 0.55 set xlabel \"MLT\" set ylabel \"Latitude\" set x2label \"Latitude\" set xtics ('12' 12, '15' 15, '18' 18, '21' 21, '24' 24, '3' 27, '6' 30, '9' 33, '12' 36) # splot 'f_aveFe_245.gnu' u 1:2:3 with pm3d, 'cont2.dat' u 1:2:3 w l splot 'tmp.1' u 1:2:3 with pm3d "; } s/;/\n/gs; s/\n\s+/\n/gs; s/^#.*$//gs; $t=$_; # print "\n----------\n$t\n-------------\n"; ssy("/bin/rm -vf tmp.gnu tmp.ps tmp.pdf"); dumpf($t,"tmp.gnu"); ssy(" gnuplot < tmp.gnu > tmp.ps"); ssy(" ps2pdf tmp.ps"); ssy("pdfcrop --margins 4 tmp.pdf; /bin/mv tmp-crop.pdf tmp.pdf"); $out=parsev($v,'out=showpdf'); if($out eq 'showpdf'){ ssy(" open tmp.pdf"); exit;} #..... maybe better w/ghostscript. convert does crappy job of rendering ps/pdf else{ ssy("convert -density 600 tmp.pdf tmp.ppm; convert -resize 720x720 -quality 100 tmp.ppm $out"); } } # # http://gnuplot-tricks.blogspot.com/2009/07/maps-contour-plots-with-labels.html # http://www.gnuplotting.org/tag/sphere/ # https://stackoverflow.com/questions/44055951/gnuplot-contour-and-heatmap-from-two-different-files-in-the-same-splot # https://stackoverflow.com/questions/18792461/gnuplot-2d-polar-plot-with-heatmap-from-3d-dataset-possible # https://stackoverflow.com/questions/18878163/gnuplot-contour-line-color-set-style-line-and-set-linetype-not-working # https://stackoverflow.com/questions/18878163/gnuplot-contour-line-color-set-style-line-and-set-linetype-not-working #------------------------------------------------- sub make_movie{ local($in,$out)=@_; #------------------------------------------------- local($a,$b,$k,$kk); $k=0; ssy("/bin/rm -vf $out"); foreach $a (glob($in)){ $kk=sprintf("%6.6d",$k++); ssy("convert -quality 100 $a tmpm_$kk.jpg"); } #.... see https://evilsoup.wordpress.com/2013/02/10/general-ffmpeg-encoding-guide-2/ # makes small but good quality mp4 ssy("ffmpeg -i tmpm_%06d.jpg -c:v libx264 -crf 20 -preset fast $out"); ssy("/bin/rm -v tmpm_??????.jpg"); } #------------------------------------------------- sub parsev{ local($v,$w)=@_; #------------------------------------------------- local($a,$b,$c,$d); ($a,$b)=split('=',$w); # print STDERR "v:|$v|\n $w ==> $a $b \n"; foreach (split(' ',$v)){ ($c,$d)=split('=',$_); # print STDERR " |$a| |$b| |$c| |$d| \n"; if($a eq $c){ return($d); } } return($b); } #------------------------------------------------- sub remote{ local($v)=@_; #------------------------------------------------- # run myself remotely local($x,$t,$rdir); $x=shift(@ARGV); if($x eq 'isremote'){ return; } $t=parsev($v,'target=none'); $rdir=parsev($v,'rdir=tmp.remoterun'); ssy("ssh $t \"mkdir -p $rdir\""); ssy("rsync -auvz --progress --size-only . ${t}:$rdir"); # ssy("ssh $t \"cd $rdir; ./$0 isremote\" >log.remote 2>&1 $z2){$z2=$rr;} } print STDERR "got all for $variable nr=$nr $z1 $z2\n"; $h=~s/\s+$//; $h=~s/\n/ /gs; @r=(@r,$h); if($vget eq $variable){ return(@r); } } } #............. collection of simple funcs in perl # now merged with plot stuff sub ssy{ local($c)=@_; print STDERR "ssy:|$c|\n"; system($c); return; } sub min{ local($a,$b)=@_;local($c); $c=$a; if($b<$a){$c=$b;} return($c); } sub max{ local($a,$b)=@_;local($c); $c=$a; if($b>$a){$c=$b;} return($c); } sub abs{ local($a)=@_;local($c); $c=$a; if($a<0){$c=-$a;} return($c); } sub log10{local($a)=@_;$y10=1.0/(log(10.0));$a=$y10*log(max($a,1.0e-15)); return($a); } sub log2{local($a)=@_;$y10=1.0/(log(2.0));$a=$y10*log(max($a,1.0e-15)); return($a); } sub basename{ local($k)=@_; local(@a)=split('\/',$k); return($a[$#a]); } sub suck{local($f)=@_; local($t,@t);$t=`cat $f`;@t=split('\n',$t); print STDERR "suck: |$f| |$#t|\n"; return(@t); } sub dumpf{ local($c,$f)=@_; open(DF,">$f"); print DF $c; close(DF); } sub dumpa{ local($f)=shift(@_); open(FF,">$f");print FF join("\n",@_); close(FF); } sub mes{ local($x)=@_; print STDERR ">$x<\n"; } sub err{ local($x)=@_; print STDERR "error $0 message: |$x|\n"; exit; } sub minmax { local(*x,*n,*mi,*ma)=@_; local($l); $mi=1.0e33; $ma=-$mi; for($l=1;$l<=$n;$l++){ if($ma<$x[$l]){$ma=$x[$l];} if($mi>$x[$l]){$mi=$x[$l];} } } sub exfile{ local($me,$ff)=@_; local($i); open(III,"<$me"); open(OOO,">$ff"); local($i)=0; while(){@a=split(); if($a[0] eq 'ENDFILE'){$i=0}; if($i==1){print OOO "$_";} if(($a[0] eq 'BEGINFILE')&&($a[1] eq $ff)){$i=1}; } close(III); close(OOO); return; } sub print_table{ local($t,$s)=@_; local($x,$y,$l,$i,$nc,@cl); for($i=0;$i<100;$i++){$cl[$i]=0; } $nc=0; foreach $l (split('\n',$t)){ $i=0; foreach $x (split($s,$l)){ $y=length($x); if($y>$cl[$i]){$cl[$i]=$y;} $i++; if($i>$nc){$nc=$i;} }} foreach $l (split('\n',$t)){ $i=0; foreach $x (split($s,$l)){ $y=$cl[$i]; $i++; $f="%-${y}s"; printf $f,$x; print STDERR " ";} print STDERR "\n"; } } sub hsvrgb { local(*h,*s,*v,*r,*g,*b)=@_; local($i,$hh,$k,$f,$m,$n); $hh=6.*$h; $i=int($hh); $f=$hh-$i; $m=$v*(1.-$s); $n=$v*(1.-$s*$f); $k=$v*(1.-$s*(1.-$f)); if($i==0){ $r=$v; $g=$k; $b=$m; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex);} if($i==1){ $r=$n; $g=$v; $b=$m; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex);} if($i==2){ $r=$m; $g=$v; $b=$k; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex);} if($i==3){ $r=$m; $g=$n; $b=$v; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex);} if($i==4){ $r=$k; $g=$m; $b=$v; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex);} $r=$v; $g=$m; $b=$n; $hex=sprintf('#%2.2x%2.2x%2.2x',256*$r,256*$g,256*$b);return($hex); } sub waitfor{ local($p,$n,$w)=@_; # name max wait/sec local($t,$m); while(1){ $t=`ps auxw`;$m=0; foreach (split('\n',$t)){if(/$p/){$m++;}} print STDERR "waitfor: $m/$n ... "; if($m<$n){ print STDERR " ..... end\n"; return; } if( -e "waitfor"){ $x=`cat waitfor`; ($n,$w)=split(' ',$x); if($w<2){$w=2;} print STDERR " new n=$n w=$w ";} print STDERR " waiting $w secs\n"; sleep($w); } } #..... tile several pdf files into one sub pdf_join{ ($a,$oo)=@_; $c=""; foreach(split(' ',$a)){ ($k,$l)=split('\:',$_); ssy("amak.pdf-stamp '($l)' 490 450 Helvetica-Bold 20 tmp$k.pdf"); ssy("pdfcrop --margins -4 tmp-99.pdf; /bin/mv tmp-99-crop.pdf tmp2-$k.pdf"); $c="$c tmp2-$k.pdf"; } ssy("pdfjam $c --nup 2x3 --outfile tmp99.pdf; pdfcrop tmp99.pdf; /bin/mv tmp99-crop.pdf $oo"); } #..... simple interface to plot 2d contours from ascii data files sub make_2dp{ local($fin,$sel,$fmima,$log,$cmap,$lab,$out)=@_; $c="amak.2dpl-20160320 fin=$fin sel=$sel col_ix=14 col_iy=15 col_vv=13 col_xx=11 col_yy=12 xybox=4:16:4:16 fmima=$fmima outlier=none cmap=$cmap cbar=yes log=$log tout=tmp.ps ch=0.55 fnt=Helvetica xax=longitude yax=latitude yaxoff=3.0 topl=~ topc=$lab topr=~ con=none ytics=none outlier=none xymax=0:999999:0:999999 xgrid=none ygrid=none mhdgrid=none points=none "; $c=~s/\n/ /gs; ssy($c); ssy("ps2pdf tmp.ps; pdfcrop --margins -3 tmp.pdf; /bin/mv -v tmp-crop.pdf $out.pdf");} #--------------------------------------------------------------------------------------- #----- time epoch routines -------------------------------------------------------- #--------------------------------------------------------------------------------------- sub epo2i{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $t=epo1($iy,$mo,32,0,0,0); $id=32;while($t>=$e){$id--;$t=epo1($iy,$mo,$id,0,0,0);} $t=epo1($iy,$mo,$id,0,0,0); $t=$e-$t; $ih=int($t/3600.0); $t=$t-(3600*$ih); $mi=int($t/60.0); $se=$t-(60*$mi); $se1=int($se); $se=$se-$se1; $se2=int((1000*$se)+0.1); $r=sprintf("%4.4d:%2.2d:%2.2d:%2.2d:%2.2d:%2.2d.%3.3d", $iy,$mo,$id,$ih,$mi,$se1,$se2); } sub epo2i_se{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $t=epo1($iy,$mo,32,0,0,0); $id=32;while($t>=$e){$id--;$t=epo1($iy,$mo,$id,0,0,0);} $t=epo1($iy,$mo,$id,0,0,0); $t=$e-$t; $ih=int($t/3600.0); $t=$t-(3600*$ih); $mi=int($t/60.0); $se=$t-(60*$mi); $r=$se; } sub epo2i_mi{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $t=epo1($iy,$mo,32,0,0,0); $id=32;while($t>=$e){$id--;$t=epo1($iy,$mo,$id,0,0,0);} $t=epo1($iy,$mo,$id,0,0,0); $t=$e-$t; $ih=int($t/3600.0); $t=$t-(3600*$ih); $mi=int($t/60.0); $r=$mi; } sub epo2i_ho{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $t=epo1($iy,$mo,32,0,0,0); $id=32;while($t>=$e){$id--;$t=epo1($iy,$mo,$id,0,0,0);} $t=epo1($iy,$mo,$id,0,0,0); $t=$e-$t; $ih=int($t/3600.0); $r=$ih; } sub epo2i_id{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $t=epo1($iy,$mo,32,0,0,0); $id=32;while($t>=$e){$id--;$t=epo1($iy,$mo,$id,0,0,0);} $r=$id; } sub epo2i_mo{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $t=epo1($iy,13,1,0,0,0); $mo=13;while($t>=$e){$mo--;$t=epo1($iy,$mo,1,0,0,0);} $r=$mo; } sub epo2i_iy{ local($e)=@_; local($t,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $t=epo1(2020,1,1,0,0,0); $iy=2020;while($t>=$e){$iy--;$t=epo1($iy,1,1,0,0,0);} $r=$iy; } sub epo2{local($e)=@_; local($iy,$mo,$id,$ih,$mi,$se)=split(':',$e); local($rr); $rr=epo1($iy,$mo,$id,$ih,$mi,$se); } sub epo1{ local($iy,$mo,$id,$ih,$mi,$se)=@_; local($t,$i1,$i2,$i3,$i4,$d1,$rr); if($iy<200){$iy+=1900;} $i1=367*$iy; $i2=int(($mo+9)/12); $i2=($i2+$iy)*7; $i2=int($i2/4); $i3=int((275*$mo)/9); $i4=100*$iy+$mo; $d1=(1.0); if(($i4-190002.5)<0.0){$d1=(-1.0);} $rr = ( $i1 - $i2 + $i3 + $id + 1721013.5 -0.5*$d1 +0.5 )*86400.0; # in seconds now $rr = $rr + 3600.0*$ih + 60.0*$mi +$se; } sub epo1t{ local($tt)=@_; local($iy,$mo,$id,$ih,$mi,$se)=split('\:',$tt); local($t,$i1,$i2,$i3,$i4,$d1,$rr); if($iy<200){$iy+=1900;} $i1=367*$iy; $i2=int(($mo+9)/12); $i2=($i2+$iy)*7; $i2=int($i2/4); $i3=int((275*$mo)/9); $i4=100*$iy+$mo; $d1=(1.0); if(($i4-190002.5)<0.0){$d1=(-1.0);} $rr = ( $i1 - $i2 + $i3 + $id + 1721013.5 -0.5*$d1 +0.5 )*86400.0; # in seconds now $rr = $rr + 3600.0*$ih + 60.0*$mi +$se; } #--------------------------------------------------------------------------------- # plotit: actually make the plot #--------------------------------------------------------------------------------- sub plotit{ local($v)=@_; local(%FMI,%FMA); %FMI=(); %FMA=(); $v=~s/\;/\n/gs; @v=split('\n',$v); $vv=$v; $vv=~s/\n/ /gs; $verb=pvv('verbose',0); print STDERR "plotit: \n$v\n\n"; $np=0;foreach(@v){if(/^\s*panel/){$np++;}} $xxa=pvv('xxa',2); $xxb=pvv('xxb',14); $yya=pvv('yya',2); $yyb=pvv('yyb',24); $T1=pvv('T1','1966:01:01:00:00:0.0'); $t1=epo1t($T1); $T2=pvv('T2','2020:01:01:00:00:0.0'); $t2=epo1t($T2); $lwd=pvv('lwd',0.05); $ch=pvv('ch',0.3); $fnt="Helvetica"; if($np<1){print STDERR "error np=0\n";exit;} $dy=($yyb-$yya)/$np; #..... loop over panels, find fmin/fmax foreach $pv (@v){ if(($pv=~/panel/)||($pv=~/addp/)){ $pn=ppv('pname','REQ');($m,$z1,$z2)=gtts(); if(exists($FMI{$pn})){$FMI{$pn}=min($z1,$FMI{$pn});}else{$FMI{$pn}=$z1;} if(exists($FMA{$pn})){$FMA{$pn}=max($z2,$FMA{$pn});}else{$FMA{$pn}=$z2;} $x=ppv('fmin','XXX'); if($x ne 'XXX'){ $FMI{$pn}=$x; } $x=ppv('fmax','XXX'); if($x ne 'XXX'){ $FMA{$pn}=$x; } print STDERR "panel |$pn| m=|$m| |$z1| |$z2|\n"; } } foreach $pn (sort(keys(%FMI))){print STDERR "mima: |$pn| |$FMI{$pn}| |$FMA{$pn}|\n";} #..... now plot panels $kp=0; plopen(); foreach $pv (@v){ if(($pv=~/panel/)||($pv=~/addp/)){ $pn=ppv('pname','REQ'); if($pv=~/panel/){ $yy=$yya+$dy*$kp;window($xxa,$xxb,$yy,$yy+$dy,$t1,$t2,$FMI{$pn},$FMA{$pn}); $x=ppv('boxlwd',0.01); $p->setlinewidth($x); yax1(ticdi($FMI{$pn},$FMA{$pn},10),0.1,1); wbox();$kp++; } print STDERR "plotting line >>|$pv|\n"; ($m,$z1,$z2)=gtts(); plcrv($FMI{$pn},$FMA{$pn}); #..... plot the line }} #...... time axis $tic=pvv('ttic',0.12); window($xxa,$xxb,$yya,$yyb,$t1,$t2,0,1);tax2($t1,$t2,$tic); pclose(); } #--------------------------------------------------------------------------------- sub fdot{ local($x,$y,$r,$c)=@_; setc($c); $p->circle({filled=>1},tx($x),ty($y),$r);} #--------------------------------------------------------------------------------- #--------------------------------------------------------------------------------- # plcv: plot a curve, lines, points, ... #--------------------------------------------------------------------------------- sub plcrv{ local($y1,$y2)=@_; local($i,@a,$x,$r,$s,$y,$fi); $x=ppv('lwd',0.012); $p->setlinewidth($x); print STDERR "plrcv setlinewidth --> $x\n"; $c=ppv('col','#000000'); setc($c); $r=ppv('rad',0.05); $s=ppv('sty','line');$k=0;@a=();$m=$#TT; print STDERR "plcrv: $x $r $s |$c|\n"; if($s eq 'line'){for($i=0;$i<=$m;$i++){$a[$k++]=tx($TT[$i]);$a[$k++]=ty($VV[$i]);} $p->polygon({filled => 0},@a); } $fi=0; if($s eq 'fdots'){ $fi=1; } if($s=~/[of]dots/){ for($i=0;$i<=$m;$i++){ $y=$VV[$i]; if(($y>=$y1)&&($y<=$y2)){ $p->circle({filled=>$fi},tx($TT[$i]),ty($y),$r);}}} black(); } #--------------------------------------------------------------------------------- # gtts: get a time series either frm file or from stored array #--------------------------------------------------------------------------------- sub gtts{ local(@a,@b,$tt,$xx,$m,$z1,$z2); $z1=1e33;$z2=-$z1; $vn=ppv('vname','REQ'); print STDERR "gtts start vn=|$vn|"; if(exists($Pa{$vn})){@a=split('\|',$Pa{$vn}); if($verb>5){print STDERR "from Pa |$#a| |$vn|" . substr($Pa{$vn},0,55) . "\n";}} else{ @a=suck($vn); } @b=split(' ',$a[0]); $istt=0; if($b[7]>1000){ $istt=1; } $m=0; foreach(@a){ @b=split(); if($istt>0){ $tt=$b[7]; } else{ $tt=epo1(@b); } if(($tt>=$t1)&&($tt<=$t2)){ $xx=$b[6]; if($xx<$z1){$z1=$xx;} if($xx>$z2){$z2=$xx;} $TT[$m]=$tt; $VV[$m++]=$xx; }} if($m<3){ print STDERR "error gtts not read enough |$m| |$vn|\n"; exit; } return( ($m,$z1,$z2) ); } #--------------------------------------------------------------------------------- # parse from variable pv #--------------------------------------------------------------------------------- sub pvv{ local($s,$d)=@_;local($s1,$s2);foreach(split(' ',$vv)){ ($s1,$s2)=split('=',$_); if($verb>5){print STDERR "pvv |$vv| |$s| |$s1| |$s2| \n";} if($s1 eq $s){return($s2);} } if($d eq 'REQ'){print STDERR "pvv cant find |$s|\n";exit;}else{return($d);}} sub ppv{ local($s,$d)=@_; local($k,$s1,$s2); if($verb>5){print STDERR "ppv start\n";} foreach $k (split(' ',$pv)){ ($s1,$s2)=split('=',$k); if($verb>5){print STDERR "ppv |$vv| |$s| |$s1| |$s2| \n";} if($s1 eq $s){return($s2);}} if($d eq 'REQ'){print STDERR "cant find |$s|\n"; exit;} else{return($d);} } #--------------------------------------------------------------------------------- # plopen start new plot #--------------------------------------------------------------------------------- sub plopen{ local($a1); $a1=pvv('paper','Letter'); $p=new PostScript::Simple(papersize=>$a1,colour=>1,eps=>1,clip=>1, direction=>"RightUp",coordorigin=>"LeftBottom",units=>"cm"); } #--------------------------------------------------------------------------------- # pclose close plot #--------------------------------------------------------------------------------- sub pclose{ local($a)=@_; $vv="$vv $a"; print STDERR "pclose pv=|$pv|"; local($of); ssy("/bin/rm -f tmppl.ps"); $of=pvv('outfile','tmp.pdf'); $p->output('tmppl.ps'); ssy("pstopdf tmppl.ps"); $a=pvv('crop','XXX'); print STDERR "crop:$a\n"; if($a ne 'XXX'){ ssy("pdfcrop --margin $a tmppl.pdf; /bin/cp tmppl-crop.pdf $of"); } else{ ssy("/bin/cp -v tmppl.pdf $of"); } $a=pvv('showit','no'); if( ($a ne 'no') && ($a ne '0') ){ ssy("open $of"); } } #--------------------------------------------------------------------------------- # crv plot a curve in window #--------------------------------------------------------------------------------- sub crv{ local($n,$z1,$z2,$li,$co,$wd)=@_; local($lws,$i,@a,$x,$y); $lws=$lwd; $lwd=$wd; $p->setlinewidth($lwd); setc($co); if($li==1){ @a=(); #.... simple line for($i=0;$i<$n;$i++){$x=tx($TT[$i]);$y=ty(max($z1,min($z2,$VV[$i]))); push(@a,$x,$y); #print STDERR "crv: $t1 $TT[$i] $t2 $x $y\n"; } $p->line($a[0],$a[1],$a[2],$a[3]); for($i=4;$i<2*$n;$i=$i+2){ $p->linextend($a[$i],$a[$i+1]); } } black(); $lwd=$lws; } #--------------------------------------------------------------------------------- # gtq quick read of time series in standard format #--------------------------------------------------------------------------------- sub gtq{ local($ff,$t1,$t2,$sc)=@_; local(@r,@a,@b,$n); @a=suck($ff); $n=0; foreach(@a){@b=split(); $t=epo1(@b); if(($t>=$t1)&&($t<=$t2)){ $TT[$n]=$t; $bb=$sc*$b[6]; $VV[$n]=$bb; $z1=min($z1,$bb); $z2=max($z2,$bb); $n++; # print STDERR " $b[6] $z1 $z2\n"; }} @r=($n,$z1,$z2); print STDERR "gtq $a[0] $a[$#a] $n $z1 $z2 $t1 $t2 sc=$sc\n"; return(@r); } #--------------------------------------------------------------------------------- # window define mapping into plot window #--------------------------------------------------------------------------------- sub window{ ($wwx1,$wwx2,$wwy1,$wwy2,$ffx1,$ffx2,$ffy1,$ffy2)=@_; print STDERR "window: ($wwx1,$wwx2,$wwy1,$wwy2,$ffx1,$ffx2,$ffy1,$ffy2) \n"; $wfx=($wwx2-$wwx1)/($ffx2-$ffx1);$wfy=($wwy2-$wwy1)/($ffy2-$ffy1);} #--------------------------------------------------------------------------------- # wbox box around plot window #--------------------------------------------------------------------------------- sub wbox{ local($lwd)=@_; $p->setlinewidth($lwd); black(); $p->box($wwx1,$wwy1,$wwx2,$wwy2); } #--------------------------------------------------------------------------------- # lcw set line color and width #--------------------------------------------------------------------------------- sub slwd{ local($a)=@_; $p->setlinewidth($lwd*$a); } sub lcw{ local($r,$g,$b,$a1); $a1=pv('lwd','XXX'); if($a1 ne 'XXX'){ $lwd=$a1; } $a1=pv('col','XXX'); print STDERR "lcw:col: |$a1|\n"; if($a1 ne 'XXX'){ ($r,$g,$b)=split(':',$a1); setc($r,$g,$b); } } #--------------------------------------------------------------------------------- # set various colors #--------------------------------------------------------------------------------- sub black{ setc(0,0,0); } sub white{ setc(1,1,1); } sub red{ setc(1,0,0);} sub darkred{ setc(0.75,0,0);} sub green{setc(0,1,0);} sub darkgreen{ setc(0,0.75,0);} sub blue {setc(0,0,1);} sub darkblue{setc(0,0,0.75);} sub aqua1{setc(0.399843,0.803606,0.666405);} sub darkorange{ setc(0.6,0.6,0.15);} sub grey1{ setc(0.8,0.8,0.9); } #--------------------------------------------------------------------------------- # setc set color either #rrggbb or R:G:B #--------------------------------------------------------------------------------- sub setc{ local($r,$g,$b)=@_; local($f,$rr,$gg,$bb,$rx,$bx,$gx); $f='3'; # print STDERR "setc in |$r|\n"; if($r=~/\#(\w{6})$/){ $g=$1; #..... hex color ($b)=sscanf("%x",$g); $rr=$b; $bx=$b%256; $b=$b/256; $gx=$b%256; $rx=int($b/256); if(0){print STDERR "setc # |$r| --> ($rx,$gx,$bx)\n";} $f='#'; } if($r=~/:/){ $f='c'; ($rx,$gx,$bx)=split(':',$r); } #..... colon separated fractions if($f eq '3'){ $rx=$r; $gx=$g; $bx=$b; } #..... otherwise 3 fractions if($f ne '#'){ $rx=int(max(0,min(255,256*$rx))); $gx=int(max(0,min(255,256*$gx)));$bx=int(max(0,min(255,256*$bx))); } if(0){print STDERR "setc: |$r| |$g| |$b| |$rr| |$gg| |$bb| |$rx| |$gx| |$bx| |$f|\n"; } $p->setcolour($rx,$gx,$bx); } #--------------------------------------------------------------------------------- # pv grab parameter from @l array #--------------------------------------------------------------------------------- sub pv{ local($a,$b)=@_; foreach(@l){($c,$d)=split('\=',$_); if($c eq $a){ return($d);}} return($b);} #--------------------------------------------------------------------------------- # ppe grab parameter from $a string with eval #--------------------------------------------------------------------------------- sub ppe{local($x,$a)=@_;local($r,$s,$w,$e); $_=" $x "; foreach $w (split(' ',$a)){ if(/\s+${w}=(\S+)\s+/){ $s=$1; $e="\$$w='$s';"; $r=eval($e); print STDERR "ppe set: |$e| |$w|-->|$r|\n"; }} } #--------------------------------------------------------------------------------- # tax2 make time axis #--------------------------------------------------------------------------------- sub tax2{ local($t1,$t2,$tic)=@_; local($x,$y,$ch0,$DT,$lev,$r,$e,@D,$h0,$h1,$h2,$d,$d0,$i,$tx,$v,$t,$m,$n); $_="1 2 3 5 10 15 20 30 60 2*60 3*60 5*60 10*60 15*60 20*60 30*60 3600 2*3600 4*3600 8*3600 12*3600 24*3600 2*24*3600 3*24*3600 4*24*3600 7*24*3600 14*24*3600 28*24*3600"; $DT=$t2-$t1; foreach $x (split){ $e=eval("\$r=$x;"); push(@D,$e); } %L=();$lev=0;$h0='';$h1='';$h2=''; foreach $d (@D){ $n=int($DT/$d); #print STDERR "dt=$DT d=$d n=$n lev=$lev d0=$d0 h0=$h0 h1=$h1 h2=$h2 \n"; if(($lev==2)&&($n<12)){ $h2=tax2a($h1,$t1,$t2,$d,$lev,$tic*1.4); } if(($lev==1)&&($n<22)){ if($d0*int($d/$d0)==$d){$h1=tax2a($h0,$t1,$t2,$d,$lev,$tic*1.2); $lev++;}} if(($lev==0)&&($n<60)){ $d0=$d; $h0=tax2a($hm,$t1,$t2,$d,0,$tic); $lev++; } } print STDERR "tic=$tic "; foreach $i (sort(keys(%L))){ $v=$L{$i}{"v"}; $t=$L{$i}{"t"}; $m=$L{$i}{"m"}; print STDERR "$v$m($t) "; mt($i,$ffy1); dd(0,$t); mt($i,$ffy2); dd(0,-$t); $tx="$v$m"; $tx=~s/60m/00m/; $tx=~s/24h/0h/; if(($tic>0)&&($m ne '')){ mt($i,$ffy1); md($adjustx,-0.5*$t-$ch+$adjusty); $ch0=$ch; $ch=$t*$ch0/(1.3*$tic);txt("$v","cc",0); $d=length("$v"); md($d*0.5*$ch,0.4*$ch); $ch=0.80*$ch; txt("$m","cc",0); $ch=$ch0; } } if($tic>0){ $h1=isoT($t1); $h2=isoT($t2); $x=0.5*($wwx1+$wwx2); $y=$wwy1-3*$ch; mw($x,$y); txt("$h1 to $h2","cc",0); } } #--------------------------------------------------------------------------------- # isoT format time string to ISO spec #--------------------------------------------------------------------------------- sub isoT{ local($t)=@_; local($iy,$mo,$id,$ih,$mi,$se,$r); $_=epo2i($t);($iy,$mo,$id,$ih,$mi,$se)=split(':'); $r=sprintf("%4.4d-%2.2d-%2.2dT%2.2d:%2.2d:%2.2d",$iy,$mo,$id,$ih,$mi,$se); return($r); } #--------------------------------------------------------------------------------- # tax2a internal time axis #--------------------------------------------------------------------------------- sub tax2a{ local($h,$t1,$t2,$dt,$lev,$tic)=@_; local($hm,$v,$i1,$i,$i2); $i1=$dt*(int($t1/$dt)-1); $i2=$dt*(int($t2/$dt)+1); for($i=$i1;$i<=$i2;$i=$i+$dt){ if(($i>=$t1)&&($i<=$t2)){ $hm='d'; $v=epo2i_id($i); if($dt<82000){ $hm='h'; $v=epo2i_ho($i); } if($dt<3600){ $hm='m'; $v=epo2i_mi($i); } if($dt<60){ $hm='s'; $v=epo2i_se($i); } #print STDERR "tax2a: i=$i dt=$dt h=$h hm=$hm v=$v\n"; if(($lev>1)&&($hm eq $h)){ return($hm); } $L{$i}{"v"}=$v;$L{$i}{"t"}=$tic;if($lev==0){$L{$i}{"m"}='';} else{ $L{$i}{"m"}=$hm;}}} return($hm);} #--------------------------------------------------------------------------------- # xax1 simple x-axis #--------------------------------------------------------------------------------- sub xax1{ local($s,$t,$fs)=@_; local($i,$i1,$i2); $i1=$s*(int($ffx1/$s)-5); $i2=$s*(int($ffx2/$s)+5); print STDERR "%% xaxis $i1 to $i2 in $s \n"; for($i=$i1;$i<=$i2;$i+=$s){ if(($i>=$ffx1)&&($i<=$ffx2)) { # print STDERR "xaxis:tic $i $ffx1 $ffx2 $s $i1 $i2\n"; mt($i,$ffy1); dd(0,$t); mt($i,$ffy2); dd(0,-$t); if($fs>=0){ mt($i,$ffy1); md($adjustx,-0.5*$t-$ch+$adjusty); txt("$i","cc",$fs); } } } } #--------------------------------------------------------------------------------- # yaxr: y-axis on right side #--------------------------------------------------------------------------------- sub yaxr{ local($s,$t,$fs)=@_; local($nn,$itxt,$i,$i1,$i2); $nn=0; mt($ffx2,$ffy1); dt($ffx2,$ffy2); # draw line again, may be colored if($s<=0){$s=1} $i1=$s*(int($ffy1/$s)-5); $i2=$s*(int($ffy2/$s)+5); for($i=$i1;$i<=$i2;$i+=$s){ if(($i>=$ffy1)&&($i<=$ffy2)) { #print STDERR "%yaxis:tic i=$i s=$s i1=$i1 i2=$i2\n"; mt($ffx2,$i); dd(-$t,0); if($fs>=0){ mt($ffx2,$i); md($ch+$adjustx,-0.5*$ch+$adjusty); if(&abs($i)<1.0e-10){$i=0;} $itxt="$i"; $itxt=~s/\0*00001$//; txt("$itxt","lc",$fs); } } ++$nn; if($nn>100){ return; } } } #--------------------------------------------------------------------------------- # y-axis log scale, fs determines density of labels, $t is tic length, $s unused #--------------------------------------------------------------------------------- sub yax1log{ local($s,$t,$fs)=@_; local($yt,$yl,$kt1,$kt2,$dol,$kp); $dol=1; $kt1=0; $kt2=0; $kp=0; foreach $w (qw/1 3 2 5 4 6 7 8 9/) { if($kt1<20*$fs){ foreach $ex (-10..10){ $yt=(10**$ex)*$w; $yl=log10($yt); if(($yl>$ffy1)&&($yl<$ffy2)){ print STDERR "yaxlog: $fs $kt1 $kt2 $dol $ffy1 $yt $yl $ffy2\n"; $kt1++; mt($ffx1,$yl); dd($t,0); mt($ffx2,$yl); dd(-$t,0); if(($fs>0)&&($dol>0)){ mt($ffx1,$yl);md(-$ch+$adjustx,-0.1*$ch+$adjusty); if( ($w!=4)&&($w<6)){ $kt2++; if($ex==0){ $x=$w; txt("$x","rc",0); $kp++; next; } if($ex==1){ $x=10*$w; $kp++;txt("$x","rc",0); next; } if($ex==2){ $x=100*$w; $kp++;txt("$x","rc",0); next; } if($ex==-1){ $x=0.1*$w; $kp++;txt("$x","rc",0); next; } if($ex==-2){ $x=0.01*$w; $kp++;txt("$x","rc",0); next; } $kp++;txt("10",'rc',0); if($w>1){ mt($ffx1,$yl);md(-$ch+$adjustx,-0.1*$ch+$adjusty); md(-2.0*$ch,0); $kp++;txt("${w}x",'lc',0); } if($ex != 0){mt($ffx1,$yl);md(-$ch+$adjustx,-0.1*$ch+$adjusty); $ch=$ch/1.4; md(0.5*$ch,0.5*$ch); $kp++;txt("$ex","rc",0); $ch=$ch*1.4; } } } } } if($kt2>4*$fs){$dol=0;} } } if($p==0){ #..... desperation, to tics yet $kt1=($ffy2-$ffy1)/7; for($yt=$ffy1+0.5*$kt1;$yt<$ffy2;$yt=$yt+$kt1){ $y1=log10($yt); mt($ffx1,$yl); dd($t,0); mt($ffx2,$yl); dd(-$t,0); mt($ffx1,$yl);md(-$ch+$adjustx,-0.1*$ch+$adjusty); txt("$yt",'rc',0); } } } #--------------------------------------------------------------------------------- # y-axis left side #--------------------------------------------------------------------------------- sub yax1{ local($s,$t,$fs)=@_; local($nn,$itxt,$i,$i1,$i2,$ych); $nn=0; if($s<=0){$s=1} $i1=$s*(int($ffy1/$s)-5); $i2=$s*(int($ffy2/$s)+5); for($i=$i1;$i<=$i2;$i+=$s){ $ych=ty($i); if(($i>=$ffy1)&&($i<=$ffy2)&&($ych+$ch<$wwy2)) { #print STDERR "%yaxis:tic i=$i s=$s i1=$i1 i2=$i2\n"; mt($ffx1,$i); dd($t*1.3,0); mt($ffx2,$i); dd(-$t*1.3,0); if($fs>=0){ mt($ffx1,$i); md(-$ch+$adjustx,-0.5*$ch+$adjusty); if(&abs($i)<1.0e-10){$i=0;} $itxt="$i"; $itxt=~s/\0*00001$//; txt("$itxt","rr",$fs); } } ++$nn; if($nn>100){ return; } } #...... unlabelled tics $s=$s/5;if($s<=0){$s=1;} $i1=$s*(int($ffy1/$s)-5); $i2=$s*(int($ffy2/$s)+5); for($i=$i1;$i<=$i2;$i+=$s){if(($i>=$ffy1)&&($i<=$ffy2)){mt($ffx1,$i);dd($t,0);mt($ffx2,$i);dd(-$t,0);}} } #--------------------------------------------------------------------------------- # calculate tic marks #--------------------------------------------------------------------------------- sub ticdi{ local($t1,$t2,$n)=@_; local($sd,$i1,$s1,$ex); $t=$t2-$t1; for($ex=-15;$ex<=15;$ex++){ foreach $s1 ( 1, 2, 5 ){ $sd=$s1*(10.0**$ex); $i1=int($t/$sd)+2; if($i1<=$n){return($sd);} } } } #--------------------------------------------------------------------------------- # plot primitives #--------------------------------------------------------------------------------- sub tx{ local($x)=@_; return($wwx1+($x-$ffx1)*$wfx); } sub ty{ local($y)=@_; return($wwy1+($y-$ffy1)*$wfy); } sub mw{ local($x,$y)=@_; $wwx9=$x; $wwy9=$y; $wwmm=0; } sub dw{ local($x,$y)=@_; if($wwmm>3){ print STDERR "error $0 wwmm=$wwmm on draw\n"; exit; } if($wwmm==0){ $p->line($wwx9,$wwy9,$x,$y); } if($wwmm==1){ $p->linextend($x,$y); } $wwx9=$x; $wwy9=$y; $wwmm=1; } sub mt{ local($x,$y)=@_; mw(tx($x),ty($y)); return(); } sub dt{ local($x,$y)=@_; dw(tx($x),ty($y)); return(); } sub md{ local($x,$y)=@_; mw($wwx9+$x,$wwy9+$y); return(); } sub dd{ local($x,$y)=@_; dw($wwx9+$x,$wwy9+$y); return(); } sub ddt{ local($x,$y)=@_; dw($wwx9+tx($x),$wwy9+ty($y)); return(); } sub cpoly{ local($x1,$x2,$y1,$y2,$rgb,$lw,$fil)=@_; # print STDERR "cpoly: ($x1,$x2,$y1,$y2,$rgb,$lw,$fil)\n"; setc($rgb); $p->setlinewidth($lw); $p->polygon({filled => $fil}, $x1,$y1, $x2,$y1, $x2,$y2, $x1,$y2, $x1,$y1); } sub txt{ local($t,$c,$a)=@_; local($yy,$vref,$href,$s); $href=substr($c,0,1); $vref=substr($c,1,1); if($wwmm>3){ print STDERR "error $0 wwmm=$wwmm on txt $t \n"; exit; } $yy=0; if($vref eq 'c'){$yy=-$ch/3;} if($vref eq 't'){$yy=-$ch;} $s=1; if(($a>0.1)&&($a<10)){ $s=$a; $a=0; } # print STDERR "txt: href=$href vref=$vref yy=|$yy| ch=|$ch| s=|$s| txt=|$t|\n"; $p->setfont($fnt,$s*$ch*72.0/2.54); $t=~s/\~/ /g; if($href eq 'l'){ $p->text({rotate=>$a,align=>'left' }, $wwx9,$wwy9+$yy,$t); } if($href eq 'c'){ $p->text({rotate=>$a,align=>'centre'}, $wwx9,$wwy9+$yy,$t); } if($href eq 'r'){ $p->text({rotate=>$a,align=>'right' }, $wwx9,$wwy9+$yy,$t); } } sub dumpa{ local($f,@a)=@_; open(DF,">$f"); foreach(@a){print DF "$_\n";} close(DF); } #------------------------------------------------------------------ sub hexcols{ local($k,$n)=@_; local($h,$s,$v,$hex); #..... k<=n #------------------------------------------------------------------ $v=0.99; $s=0.4+0.59*($k/$n); $h=$k/$n; $h=$h%1.0; $hex=hsvh($h,$s,$v); print STDERR "hexcol: k=$k n=$n h=$h s=$s v=$v $hex\n"; return($hex); } #------------------------------------------------------------------ sub hsvh{ local($h,$s,$v)=@_; local($i,$hh,$k,$f,$m,$n,$fo); #------------------------------------------------------------------ $fo='#%2.2x%2.2x%2.2x'; $hh=6.*$h; $i=int($hh);$f=$hh-$i;$m=$v*(1.-$s);$n=$v*(1.-$s*$f);$k=$v*(1.-$s*(1.-$f)); print STDERR "h=$h m=$m n=$n k=$k\n"; if($i==0){$r=$v;$g=$k;$b=$m;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} if($i==1){$r=$n;$g=$v;$b=$m;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} if($i==2){$r=$m;$g=$v;$b=$k;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} if($i==3){$r=$m;$g=$n;$b=$v;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} if($i==4){$r=$k;$g=$m;$b=$v;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} $r=$v;$g=$m;$b=$n;$hex=sprintf($fo,256*$r,256*$g,256*$b);return($hex);} #------------------------------------------------------------------ sub epo3i{ #... faster, but does not go back as far, see https://en.wikipedia.org/wiki/Julian_day #------------------------------------------------------------------ local($js)=@_; local($J,$t,$e,$f,$g,$h,$iy,$mo,$id,$ih,$mi,$se,$se1,$se2,$r); $J=int(0.5+($js/86400)); #.... comes in seconds and we are 1/2 day off????? $t=4*$J+274277; $t=int($t/146097)*3; $t=int($t/4)-38; $f=$J+1401+$t; $e=4*$f+3; $g=($e%1461); $g=int($g/4); $h=5*$g+2; $id=($h%153); $id=int($id/5)+1; $mo=int($h/153)+2; $mo=($mo%12)+1; $iy=int($e/1461)-4716+int((14-$mo)/12); $t=$js-$J*86400+43200; $ih=int($t/3600); $t=$t-3600*$ih; $mi=int($t/60); $se=$t-(60*$mi); $se1=int($se); $se=$se-$se1; $se2=int((1000*$se)+0.1); $r=sprintf("%4.4d:%2.2d:%2.2d:%2.2d:%2.2d:%2.2d.%3.3d",$iy,$mo,$id,$ih,$mi,$se1,$se2);return($r); } #------------------------------------------------------------------ # tser_ave: average ts, optional linear interpolate over gaps #------------------------------------------------------------------ sub tser_ave{ local($t1,$t2,$ave,$fill,@a)=@_; if($#a<5){ print STDERR "error tser_ave{ local($t1,$t2,$ave,$fill not enoght data\n"; exit; } local($ist,$n,$s,$t,$i,@tt,@vv,$k,$tt1,$tt2,$s1,$s2,@vn,@tn); $ist=0; @b=split(' ',$a[0]); if($b[7]>1000){$ist=1;} $n=$#a+1; print STDERR "tser_ave ist=$ist |$a[0]| n=$n ave=$ave \n"; for($i=0;$i<$n;$i++){@b=split(' ',$a[$i]);if($ist==0){$t=epo1(@b);}else{$t=$b[7];} if(($t>=$t1)&&($t<=$t2)){ $vv[$k]=$b[6]; $tt[$k]=$t; $k++; }} $tt1=$t1; $s=int($tt1/$ave); $tt1=($s-0.5)*$ave; $i=0; $n=0; @tn=(); @vn=(); for($tt1=($s-0.5)*$ave;$tt1<$t2;$tt1=$tt1+$ave){ $tt2=$tt1+$ave; $tn[$n]=0.5*($tt1+$tt2); $vn[$n]=1e33; $s1=0;$s2=0;while($i<$k){$t=$tt[$i]; if($t>$tt2){last;} $s1=$s1+$vv[$i]; $s2++; $i++; } if($s2>0){ $vn[$n]=$s1/$s2; } $n++; } print STDERR "finish ave n=$n \n"; # debug @tt=();for($i=0;$i<$n;$i++){$_=epo3i($tn[$i]);s/:/ /g;$t="$_ $vn[$i] $tn[$i]";push(@tt,$t);} return(@tt); if($fill==0){@tt=();for($i=0;$i<$n;$i++){if($vn[$i]<0.9e33){ $_=epo3i($tn[$i]);s/:/ /g;$t="$_ $vn[$i] $tn[$i]";push(@tt,$t);}} return(@tt); } #...... with fill $s1=0;$s2=0; #... go forward, remember last valid, then backwards, takes care of ends too for($i=0;$i<$n;$i++){if($vn[$i]<0.9e33){$s1=$vn[$i];$s2=0;$tt[$i]=0;}else{$s2++;$tt[$i]=$s2;$vv[$i]=$s1;}} # for($i=0;$i<$n;$i++){print STDERR "forward $i $tn[$i] $vn[$i] $tt[$i] $vv[$i]\n"; } for($i=$n-1;$i>=0;$i--){if($vn[$i]>0.9e33){ if($s2==0){$s2=$tt[$i];} $vn[$i]=$s1-($s2-$tt[$i]+1)*($s1-$vv[$i])/($s2+1);} else{$s2=0;$s1=$vn[$i]; } } # for($i=0;$i<$n;$i++){print STDERR "backward $i $tn[$i] $vn[$i] $tt[$i] $vv[$i]\n"; } #..... create new array @tt=();for($i=0;$i<$n;$i++){$_=epo3i($tn[$i]);s/:/ /g;$t="$_ $vn[$i] $tn[$i]";push(@tt,$t);} return(@tt); } 1;