sub calc_rez {
    kinematics();
    $L_prime = $L*cos($lambda);
    $dk_res = $epsilon/$L_prime**2*sqrt(720/($n_m + 4));
    $pt = $p*cos($lambda);
    $dk_ms = 0.016/($L*$p*$beta*cos($lambda)**2)*sqrt($n_rl);
    $dk = sqrt($dk_res**2 + $dk_ms**2);
    $k = 0.3*$B/$pt;
    $dp_over_p = $dk/$k;
    print OUT "p $p dp/p $dp_over_p\n";
    #print "p $p dk_ms $dk_ms dk_res $dk_res dk $dk pt $pt k $k dp/p $dp_over_p\n";
}
sub calc_angle {
    kinematics();
    $theta0 = 0.0136/($beta*$c*$p)*sqrt($n_rl + $n_rl_front)*(1.0 + 0.038*log($n_rl + $n_rl_front));
    $dtheta_ms = $theta0/sqrt(3.0);
    $dtheta_res = $epsilon/$L*sqrt(12*($n_m - 1)/($n_m*($n_m + 1)));
    $dtheta = sqrt($dtheta_res**2 + $dtheta_ms**2);
    #print "beta c p n_rl n_rl_front theta0 dtheta_ms dtheta_res\n";
    #print "$beta $c $p $n_rl $n_rl_front $theta0 $dtheta_ms $dtheta_res\n";
    print OUTA "p $p dtheta $dtheta\n";
}
sub kinematics {
    $c = 1;
    $E = sqrt($p**2 + $m**2);
    $gamma = $E/$m;
    $beta = sqrt(1.0 - 1/$gamma**2);
}
1;