#!/usr/bin/perl
$file=$ARGV[0];
my $pi = 3.14159265358979;
my @data;
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );
for my $i (0..$#data)
{
push @pri, $data[$i][0];
push @sec, $data[$i][1];
}
for my $i (0..$#data-1)
{
$dx[$i]=$data[$i][1]-$data[$i+1][1]; #col 2 diff
if($dx[$i] != 0)
{
push @ind,$i,$i+1;
}
}
push @ind1, 0,@ind, $#data;
for my $i (0..$#data)
{
if($ind1[$i] != 0)
{
push @pria, $pri[$ind1[$i]];
push @seca, $sec[$ind1[$i]]
}
}
push @x1,$pri[0],@pria;
push @y1,$sec[0],@seca;
### try this option first
#$x0=0.1+(min(@pri)+max(@pri))/2;
#$y0=0.2+(min(@sec)+max(@sec))/2;
### if you see 'saw' shape, privide center point manually
$x0=800.1;
$y0=3530.2;
#calculate azimuth of each point relative to the center point
for my $i (0..$#x1)
{
$dy[$i]=abs($y1[$i]-$y0);
$dx[$i]=abs($x1[$i]-$x0);
$theta[$i]=rad_to_deg(atan($dy[$i]/$dx[$i]));
if ($x1[$i] < $x0 && $y1[$i]>$y0 )
{
$thetacor[$i]=180-$theta[$i];
}elsif ($x1[$i] < $x0 && $y1[$i]<$y0 )
{
$thetacor[$i]=180+$theta[$i];
}elsif ($x1[$i] > $x0 && $y1[$i]<$y0 )
{
$thetacor[$i]=360-$theta[$i];
}else
{
$thetacor[$i]=$theta[$i];
}
push(@{$dataku[$i]}, $x1[$i],$y1[$i],$thetacor[$i]); #push multidimensional array
}
open OUT," > out.txt" or die "$!\n";
my @sorted = sort { $a->[2] <=> $b->[2] } @dataku; ##sort by col 2 (azimuth)
#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]}-1)
{
print OUT "$sorted[$i][$j]","\t";
}
print OUT "\n";
}
for my $j (0..$#{$sorted[0]}-1) #close polygon
{
print OUT "$sorted[0][$j]","\t";
}
print OUT "\n";
sub max {
splice(@_, ($_[0] > $_[1]) ? 1 : 0, 1);
return ($#_ == 0) ? $_[0] : max(@_);
}
sub min {
splice(@_, ($_[0] > $_[1]) ? 0 : 1, 1);
return ($#_ == 0) ? $_[0] : min(@_);
}
sub deg_to_rad { ($_[0]/180) * $pi }
sub rad_to_deg { ($_[0]/$pi) * 180 }
sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }
sub acos { atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ) }
sub tan { sin($_[0]) / cos($_[0]) }
sub atan { atan2($_[0],1) };
########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP} <<'__GNUPLOT__';
set multiplot; # get into multiplot mode
set size 1,0.5;
set origin 0.0,0.5;
set xrange[700:1150]
set yrange[3480:3570]
plot 'PRIM_SEC2.TXT' using 1:2 with points pointtype 0 pointsize 1 title ""
set origin 0.0,0.0;
set xrange[700:1150]
set yrange[3480:3570]
plot 'out.txt' using 1:2 with lines title ""
__GNUPLOT__
close $GP;
Kode di bawah ini adalah versi 2:
#!/usr/bin/perl
$file=$ARGV[0];
my @data;
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row;
}
close( FILE );
for my $i (0..$#data)
{
push @pri, $data[$i][0];
push @sec, $data[$i][1];
}
for my $i (0..$#data-1)
{
$dx[$i]=$data[$i][1]-$data[$i+1][1]; #col 2 diff
if($dx[$i] != 0)
{
push @ind,$i,$i+1;
}
}
push @ind1, 0,@ind, $#data;
for my $i (0..$#data)
{
if($ind1[$i] != 0)
{
push @pria, $pri[$ind1[$i]];
push @seca, $sec[$ind1[$i]]
}
}
push @x1,$pri[0],@pria;
push @y1,$sec[0],@seca;
for my $i (0..$#x1/2)
{
push @{$left[$i]},$x1[$i*2],$y1[$i*2]; #push multidimensional array
push @{$right[$i]},$x1[($i*2)+1],$y1[($i*2)+1]; #push multidimensional array
}
my @sortedleft= sort { $a->[1] <=> $b->[1] } @left; ##sort accending by col 2
my @sortedright = sort { $b->[1] <=> $a->[1] } @right; ##sort decending by col 2
push @sorted, @sortedleft, @sortedright;
open OUT," > out.txt" or die "$!\n";
#Print output
for my $i (0..$#sorted)
{
for my $j (0..$#{$sorted[0]})
{
print OUT "$sorted[$i][$j]","\t";
}
print OUT "\n";
}
for my $j (0..$#{$sorted[0]}) #close polygon
{
print OUT "$sorted[0][$j]","\t";
}
print OUT "\n";
########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP} <<'__GNUPLOT__';
set multiplot; # get into multiplot mode
set size 1,0.5;
set origin 0.0,0.5;
set xrange[700:1150]
set yrange[3480:3570]
plot 'PRIM_SEC.TXT' using 1:2 with points pointtype 0 pointsize 1 title ""
set origin 0.0,0.0;
set xrange[700:1150]
set yrange[3480:3570]
plot 'out.txt' using 1:2 with lines title ""
__GNUPLOT__
close $GP;
No comments:
Post a Comment