view release on metacpan or  search on metacpan

lib/Physics/Ellipsometry/VASE.pm  view on Meta::CPAN

            $delta_model -= $correction;
        }
        
        $ym .= $y_model;

        # Numerical partial derivatives via finite differences
        my $np  = $par->nelem;
        for my $i (0 .. $np - 1) {
            my $par_h = $par->copy;
            my $p_i = $par->slice("($i)")->sclr;
            my $eps = abs($p_i) * ($self->{deriv_step}) + 1e-10;
            $eps = $self->{min_deriv_step} if $eps < $self->{min_deriv_step};
            $par_h->slice("($i)") += $eps;
            my $y_pert = &$model($par_h, $x_data);
            if ($circular) {
                my $dm = $y_pert->slice("$npts:" . (2*$npts-1));
                my $diff2 = $dm - $delta_data;
                $dm -= 360.0 * rint($diff2 / 360.0);
            }
            $dyda->slice(",($i)") .= ($y_pert - $ym) / $eps;
        }
    };

    my ($ym, $finalp, $covar, $iters) = lmfit(
        $x_fit, $y_data, $sigma, $fit_func, $initial_params,
        {Maxiter => $self->{maxiter}, Eps => $self->{eps}}
    );

    $self->{covar} = $covar;
    $self->{iters} = $iters;
    $self->{ym}    = $ym;

    return $finalp;
}

# Calculate MSE (WVASE convention: sqrt(χ²/(2N-M)))
sub mse {
    my ($self, $fit_params, %opts) = @_;
    my $data  = $self->{data};
    my $model = $self->{model};
    my $nparams = $opts{nparams} // $fit_params->nelem;

    my $x_data = $data->(0:1,:)->xchg(0,1);
    my $npts   = $data->getdim(1);
    my $y_data = $data->((2),:)->flat->append($data->((3),:)->flat);

    my $y_fit = &$model($fit_params, $x_data);

    # Apply circular delta alignment for MSE calculation
    if ($self->{circular_delta}) {
        my $delta_data  = $data->((3),:)->flat;
        my $delta_model = $y_fit->slice("$npts:" . (2*$npts-1));
        my $diff = $delta_model - $delta_data;
        $delta_model -= 360.0 * rint($diff / 360.0);
    }

    my $chi2 = sum(($y_data - $y_fit)**2)->sclr;
    return sqrt($chi2 / (2*$npts - $nparams));
}

# Plot raw data with model fit overlay
sub plot {
    my ($self, $fit_params, %opts) = @_;
    require PDL::Graphics::Gnuplot;

    my $data  = $self->{data};
    my $model = $self->{model};

    my $wavelength = $data->((0),:)->flat;
    my $angles     = $data->((1),:)->flat;
    my $psi_data   = $data->((2),:)->flat;
    my $delta_data = $data->((3),:)->flat;

    # Evaluate model at fitted parameters
    my $x_data  = $data->(0:1,:)->xchg(0,1);
    my $y_model = &$model($fit_params, $x_data);
    my $npts    = $wavelength->nelem;
    my $psi_fit   = $y_model->slice("0:" . ($npts - 1));
    my $delta_fit = $y_model->slice("$npts:" . (2 * $npts - 1));

    my $output = $opts{output};
    my $title  = $opts{title} // 'VASE Fit Results';

    # Find unique angles for grouping
    my @unique_angles = sort { $a <=> $b }
                        do { my %s; grep { !$s{$_}++ } list $angles };

    # Color palette for multiple angles
    my @colors = ('#0072B2', '#D55E00', '#009E73', '#CC79A7', '#F0E442',
                  '#56B4E9', '#E69F00', '#000000');

    # Select terminal and construct gpwin
    my $gp;
    if ($output) {
        my ($term, @topts);
        if    ($output =~ /\.png$/i) { $term = "pngcairo"; @topts = (size => [900,700,"px"]) }
        elsif ($output =~ /\.pdf$/i) { $term = "pdfcairo"; @topts = (size => [7,5.5,"in"]) }
        elsif ($output =~ /\.svg$/i) { $term = "svg";      @topts = (size => [900,700,"px"]) }
        elsif ($output =~ /\.eps$/i) { $term = "epscairo" }
        else                         { $term = "pngcairo"; @topts = (size => [900,700,"px"]) }
        $gp = PDL::Graphics::Gnuplot::gpwin($term, output => $output, enhanced => 1, @topts);
    } else {
        $gp = PDL::Graphics::Gnuplot::gpwin(enhanced => 1);
    }

    # Multiplot: Psi on top, Delta on bottom
    $gp->multiplot(layout => [1, 2], title => $title);

    # Build plot curves grouped by angle
    my (@psi_curves, @delta_curves);
    for my $ai (0 .. $#unique_angles) {
        my $ang = $unique_angles[$ai];
        my $mask = ($angles == $ang);
        my $idx = which($mask);
        my $wl   = $wavelength->index($idx);
        my $psid = $psi_data->index($idx);
        my $deld = $delta_data->index($idx);
        my $psif = $psi_fit->index($idx);
        my $delf = $delta_fit->index($idx);
        my $col  = $colors[$ai % scalar @colors];
        my $label = sprintf("%.1f{/Symbol \260}", $ang);