PEYE-CAA FaultMap Generation

sub FaultMapShorts {
    my $layer = shift;
    return FaultMapLayers( \&extra_material_critical_area, [$layer], @_ );
}

sub FaultMapLayers {
    my $functionref = shift;
    my $reflayers   = shift;
    my $minspace    = shift;    # The minimum spacing
    my $plotlistref = shift;

    my $power    = 3.0;
    my $rounding = shift;
    $rounding = 3 unless ( defined $rounding );
    my @values = (
        0.55, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
        1.8,  2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0
    );                          #defect rad
    my $scalefact =
      ( $minspace * 0.6 )**$power;    #scale to make smallest defect black
    my $boundary;

    if ( defined &main::my_peye_region ) {
        $boundary = main::my_peye_region();
    }
    else {
        $boundary = region();
    }
    @fault_graphs = ();   # reset coords array -- used for graphs/tables

    #assume 1/x^power defect size distribution.
    #assume 1 defect per cm^2 of size $distance or greater
    # k/(power-1)*x^(power-1) = 1
    my $k = ( ( $power - 1 ) * $minspace**( $power - 1 ) ) /
      100000000;          # per square micron

    my $faults = 0.0;
    my $count  = 1;
    my ( $last, $val, $lastval, $lastcrit );
    foreach my $val (@values) {
        $val = $minspace * $val;
        my $crit;
        if ( $rounding == 0 ) {
            $crit = And( $boundary, $functionref->( @$reflayers, $val ) );
        }
        else {
            $crit =
              And( $boundary, $functionref->( @$reflayers, $val, $rounding ) );
        }

        # plot only every second and those that can be seen
        if ( $count <= 12 and !( $count % 2 ) ) {
            unshift( @$plotlistref, ($crit) );
            my $num = $count / 2;
            $crit->{name}  = "EYE_FM$num";    #layer names for slam/Virtuoso
            $crit->{color} = "Black";

            # 1/(defect size)**3
            $crit->{shade} = 1 - $scalefact / ($val)**$power;
        }
        $count++;
        my $critarea = area($crit);
        $val = $val * 2;                      # was radius -- now as diameter
        if ($lastval) {
            # calc straight line as y = ax + b 
            my $a = ( $critarea - $lastcrit ) / ( $val - $lastval );
            my $b = $lastcrit - $lastval * $a;

            # integrate with CA * 1/x**power
            $faults=$faults + $a/(-$power+2)*$val**(-$power+2)  
                            - $a/(-$power+2)*$lastval**(-$power+2)
                            + $b/(-$power+1)*$val**(-$power+1)  
                            - $b/(-$power+1)*$lastval**(-$power+1);
        }
        $lastval  = $val;
        $lastcrit = $critarea;
        push( @fault_graphs, # store data for graphs
            [ $val, $critarea, $critarea * $k / ( $val**($power) ) ] );
    }
    $faults = $k * $faults;

    # how many cells will fit in 1cm^2
    $faults = $faults * 100000000 / area($boundary)
      if ( area($boundary) != 0 );    # avoid float point exception
    return $faults;
}