function [f,df,ddf,x]=rsgeng(N,rL,h,lc,seed);
%RSGENG generates 1D Gaussian random rough surfaces with Gaussian Spectrum.
%
% [f,df,x]=rsgeng(N,rL,h,lc,seed)
%
% INPUT:
%
% N=total number of sample points
% rL=rough surface length
% h=rms height
% lc=correlation length
% seed=seed of random number generator
%
% OUTPUT:
%
% f=rough surface profile
% df=df/dx
% x=sample points on the surface
%
% -- Part of the Electromagnetic Wave MATLAB Library (EWML) --
% <http://www.emwave.com/>
% Original: L. Tsang, 1998.
randn('seed',seed);
y=randn(N,1);
for n=1:(N/2-1);
bh(n)=(y(2*n-1)+i*y(2*n))/sqrt(2);
end;
bhc=conj(bh);
bhf=fliplr(bhc);
bi=[bh y(N-1) bhf y(N)];
kx=2*pi*[-N/2+1:1:N/2]/rL;
y1=sqrt(wk(kx,h,lc));
y=y1*sqrt(2*pi*rL);
b=y.*bi;
xs=[b(N/2+1:1:N) b(1:1:N/2)];
xt=[xs(N),xs(1:1:N-1)];
ft=ifft(xt,N);
ft=ft*N/rL;
fs=[ft(2:1:N),ft(1)];
f=[fs(N/2+1:1:N) fs(1:1:N/2)];
f=real(f);
dx=rL/N;
x=[-N/2+1:1:N/2]*dx;
n=2:N-1;
df1=(f(n+1)-f(n-1))/(2*dx);
df=[(f(2)-f(N))/(2*dx),df1,(f(1)-f(N-1))/(2*dx)];
ddf1=(df(n+1)-df(n-1))/(2*dx);
ddf=[(df(2)-df(N))/(2*dx),ddf1,(df(1)-df(N-1))/(2*dx)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gaussian spectral density %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=wk(kx,h,lc)
y=h^2*lc*exp(-(kx*lc*0.5).^2)/(2*sqrt(pi));