complete.zip_rough surface_粗糙_粗糙面_蒙特卡罗 粗糙_随机粗糙面建模

function [R,brdf,scat,shad,w] = goa3D(f,x,y,thd,phd,n1,k1,n2,k2,nfrp,dstrb,br) % % [R,brdf,scat,shad,w] = goa3D(f,x,y,thd,phd,n1,k1,n2,k2,nfrp,dstrb,br) % % calculates the directional-hemispherical reflectance for light % incident from a medium with complex refractive index n1 + i*k1 on a 2-d % rough surface f(x,y) with index n2 + i*k2 with the angles of incidence % thd (polar) and phd (azimuth) (both in degrees) using the % geometric optics approximation (GOA) with nfrp^2 first reflection points. % Beam radius is set through br parameter while beam energy distribution is % controlled via parameter dstrb ('gauss' or 'uni'). % % Input: f - surface height % x - surface point % y - surface point % thd - (global) angle of incidence in degrees % n1 - refractive index of medium 1 % k1 - extinction coefficient of medium 1 % n2 - refractive index of medium 2 % k2 - extinction coefficient of medium 2 % nfrp - number of first reflection points (rays) along square % side % dstrb - beam energy distribution (type 'uni' for uniform and 'gauss' % for gaussian) % br - beam radius % % Output: R - directional-hemispherical reflectance % brdf - bidirectional reflection distribution function % resolved into first, second and higher order scattering % scat - mean number of scattering events per incident ray % shad - amount of shadowing (percentage) % w - percentage of rays removed due to warning signals given % during the simulation (indication of too few surface % points) % % Last updated: 2010-09-30 (David Bergstr�m) % format long nsp = length(x); % number of surface points along edge % check input if nfrp > nsp error('Too many rays! Unable to fit %g rays on surface with %g surface points.',nfrp,nsp); end; if br > (x(end) - x(1)) br = x(end) - x(1); % fit beam on surface end; % initializations theta = thd*pi/180; % Polar (zenith) AOI in radians phi = phd*pi/180; % Azimuth AOI in radians endpx = length(x); endpy = length(y); % endpoints brdf = zeros(3,91,360); % multiple scattering resolved brdf (1st order) totalenergy = 0; % total incident energy on surface reflectedenergy = 0; shp = zeros(nsp,nsp); % number of shadowed points scp = zeros(nsp,nsp); % number of scattering points per incoming ray (for every 1st refl. point) w = zeros(nsp,nsp); % number of warnings due to invalid AOI DY = (y(end)-y(1))/length(y); % sample distance slask = conv2(f(:,:),[1 0 -1],'same')/2/DY; ddy = zeros(nsp); ddy(2:end-1,2:end-1) = slask(2:end-1,2:end-1); DX = (x(end)-x(1))/length(x); % sample distance slask = conv2(f(:,:),[1 0 -1]','same')/2/DX; ddx = zeros(nsp); ddx(2:end-1,2:end-1) = slask(2:end-1,2:end-1); normals_y = -ddy; normals_x = -ddx; normals_z = ones(nsp,nsp); norm_of_normal = sqrt(normals_x.^2+normals_y.^2+normals_z.^2); % norm of surface normal normals_x = normals_x ./ norm_of_normal; % normals_y = normals_y ./ norm_of_normal; % normalizations normals_z = normals_z ./ norm_of_normal; % XX = x; YY = y; ZZ = f; normals_xx = normals_x; normals_yy = normals_y; normals_zz = normals_z; %%%%% Ray tracing algorithm starts here %%%%% for nyi = 1:nfrp; % loop over y ny = round(nsp/nfrp*nyi); % translate index to first reflection point for nxi = 1:nfrp; % loop over x nx = round(nsp/nfrp*nxi); % translate index to first reflection point % reset surface x = XX; y = YY; f = ZZ; normals_x = normals_xx; normals_y = normals_yy; normals_z = normals_zz; reflection_point = [x(nx) y(ny) f(nx,ny)]; incomingray = -[sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta)]; if shadow(x,y,f,nx,ny,theta,phi) == 0 % Test for shadowing scp(nx,ny) = 1; % scattering event found w(nx,ny) = 0; incomingenergy = incidentenergy(incomingray,normals_x,normals_y,normals_z,x,y,ddx,ddy,nx,ny,dstrb,br); % calculation of energy of incident ray incomingenergy_orig = incomingenergy; % transform incident ray and energy to scattered ray and energy [outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals_x,normals_y,normals_z,nx,ny,n1,k1,n2,k2); scattered = 0; if outgoingenergy > incomingenergy w(nx,ny) = 1; outgoingenergy = incomingenergy; outgoingray = incomingray; scp(nx,ny) = 0; scattered = 1; else totalenergy = totalenergy + incomingenergy; % add to total incident energy end; if scattered == 0 % not a weird behaving ray % check if new reflection point is available [nextreflpoint,oobpoint] = findnewpoint(x,y,f,nx,ny,outgoingray); if nextreflpoint(1) == 0 && nextreflpoint(2) == 0 % no new relection point found [scdir_theta, scdir_phi] = scatdir(outgoingray); if scdir_theta <= pi/2 % definetely scattered scattered = 1; else % out-of-bounds end; end; end; while scattered == 0 if oobpoint(1) == 0 && oobpoint(2) == 0 % ray not out-of-bounds reflection_point = [x(nextreflpoint(1)) y(nextreflpoint(2)) f(nextreflpoint(1),nextreflpoint(2))]; scp(nx,ny) = scp(nx,ny) + 1; % new scattering event found incomingenergy = outgoingenergy; incomingray = outgoingray; % transform incident ray and energy to scattered ray and % energy [outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals_x,normals_y,normals_z,nextreflpoint(1),nextreflpoint(2),n1,k1,n2,k2); % error for grazing incident rays (decrease error by increasing number of % surface points) if outgoingenergy > incomingenergy w(nx,ny) = 1; outgoingenergy = incomingenergy; outgoingray = incomingray; scp(nx,ny) = 0; totalenergy = totalenergy - incomingenergy_orig; scattered = 1; end; % check if new reflection point is available [np,oobpoint] = findnewpoint(x,y,f,nextreflpoint(1),nextreflpoint(2),outgoingray); else %bo(rp) = bo(rp) + 1; switch oobpoint(1) case endpx switch oobpoint(2) case 1 % +x -y % mirrored surface for beams out-of-bounds to lower right and upper left [np,oobpoint,ro] = findmirrorpoint(x,y,rot90(f.',2),oobpoint(1),oobpoint(2),outgoingray,8,reflection_point); % new surface f = rot90(f.',2); normals_x = -rot90(normals_x.',2); normals_y = -rot90(normals_y.',2); normals_z = rot90(normals_z.',2); case endpy % +x +y % mirro