sipangzi02.m

用monte carlo 方法计算圆周率pi,采用的是随机投点法。以及另一个应用程序

%  Spheremc.m
%  Program to Determine Fraction of Area in Latitude Bands on a Sphere
%  Tom Huber  June 25, 1996
Theta1 = 0;
Theta2 = 90;
NSubDiv = 9; % Nine Subdivisions of 10 Degrees Each
dTh = (Theta2-Theta1)/NSubDiv;  % Width of Each Division (10 Degrees)
ThLow = Theta1:dTh:Theta2-dTh;  % Lower Limit for Each Region ( 0,10,20..80)
ThHigh = Theta1+dTh:dTh:Theta2; % Upper Limit for Each Region (10,20,30..90)
Nrand = 8192;    % Number of Points for Student Edition of Matlab
Nmax  = input('How Many Loops of 8192 Values Each ');
NTrand = 0;               % Initialize Total number of Points Generated
NGoodPts = 0;             % Initialize Total number of Points on Sphere
NZone = zeros(1,NSubDiv); % Initialize Number in each zone
T0 = clock;               % Keep track of CPU time (for reference purposes)
for nloops=1:Nmax
   Xrand = rand(1,Nrand);    % Generate XYZ Points in space
   Yrand = rand(1,Nrand);
   Zrand = rand(1,Nrand);
   Rrand = Xrand.^2 + Yrand.^2 + Zrand.^2;  % Find distance from origin
   CheckValue = Rrand<=1.01 & Rrand>=.99;   % See if on surface of sphere
   NGoodPts = NGoodPts + sum(CheckValue);   % Keep track of total on surface
   Lat = asin(Zrand)*180/pi;                % Find the Latitude of Each point
   for i=1:NSubDiv                          % Sweep through all latitudes
       NZoneCheck = Lat < ThHigh(i) & Lat >= ThLow(i); % Check if in latitude
       NZoneCheck = NZoneCheck .* CheckValue;          %  and on surface
       NZone(i) = NZone(i) + sum(NZoneCheck);          % If so, add to sums
   end
   NTrand = NTrand + Nrand;                 % Total number of Points Generated
end
T0 = clock - T0;                            % CPU Time at end of program
disp(['Total Generated: ' num2str(NTrand) ' Good Pts: '...
    num2str(NGoodPts)  ' Seconds: ' num2str(T0)]);
fLatitude = NZone/NGoodPts;
fError = fLatitude./sqrt(NZone);
fActual = sin(ThHigh*pi/180.)-sin(ThLow*pi/180.);
disp('Summary for Zones');
disp('Lower Angle, Upper Angle, Simulated Fraction in Band, Uncertainty,');
disp('   Actual Fraction (using Calculus)');
disp([ ThLow' ThHigh' fLatitude' fError' fActual']);