%plankbbmesh.m
%Create a 3D printable Plot of a continuous Plank Blackbody spectra range between two
%temperatures and wavelengths.
%Steven Sahyun November 3, 2015
clc;
clf;
clear;

hc = 1240;%eV*nm Plank's constant * c.
kb = 8.617e-5; %eV/K Boltzman's constant.

%set up x vector for plot grid
xplotwidth = 100;
xmin = 0;
xmax = xmin+xplotwidth;
xstep = (xmax-xmin)/100; %100 data points
x = xmin:xstep:xmax;
x0 = 0*x+1;

%set up y vector for plot grid.
yplotwidth = 100;
ymin = 0; %Change the start point in terms of # of standard deviations
ymax = yplotwidth; %Change the end point in terms of # of standard deviations
ystep = (ymax-ymin)/100; %100 data points
y = ymin:ystep:ymax;
y0 = 0*y+1;

%create mesh grid
meshgrid(x,y);

%create a z matrix from x and y and fill with 0
z = y0'*x0;
z(:,:)=0;

%set low and high wavelengths to plot and set for 100 data points between
lambdalow = 10; %nm. can't have a wavelength of 0 nm.
lambdahigh = 1500;%nm
dlambda = (lambdahigh-lambdalow)/100;
lambda = lambdalow:dlambda:lambdahigh; %wavelength in nm

%set low and high temperatures and step range for 100 data points between 
Tlow = 4000; %low temp in K
Thigh = 7000;% high temp in K
dT = (Thigh-Tlow)/100; %100 steps
Tstepmax = Tlow+100*dT;

%This is just to get it to look nice on the plot withouth having to set 
%the axis.
%zscale = 1e17;
zscale = 60;%set height to something the printer can deal with.

%fill the z matrix with BBspectra graphs for each row changing the temp by
%dT

for Tstep = 1:101;
    T = Tlow+dT*Tstep;%go from Tlow to Thigh in 100 steps of size dT
    plankBB = 1./((lambda.^5).*(exp(hc./(kb*T*lambda))-1));
 %   Area = trapz(lambda, plankBB); %find area under curve
 %   plankBB_norm = plankBB./Area; %normalizes so area under curve = 1
    indexmax = find(max(plankBB) == plankBB);%find max value on curve
    plankBB_max = plankBB(indexmax);%scale curve so all have max value of 1
    plankBB_norm = plankBB./plankBB_max;
    A = zscale*y0'*plankBB_norm; %create a temporary matrix
  
    z(Tstep,:)=A(Tstep,:); %copy only row of interest into the z matrix
end

%plot the BB spectrum
figure(1)
plot(lambda, plankBB);

%show a 3D plot of the BB spectrum
figure(2)
surf(lambda, Tlow+y*dT, z);

%Convert to a 3D printable file; You need to have stlwrite.m and surf2solid.m by
%Sven Holcombe in the same directory as this file.
%See: http://jmumakerlab.blogspot.com/2013/11/exporting-stl-from-matlab.html
% http://www.mathworks.com/matlabcentral/fileexchange/42876-surf2solid-make-a-solid-volume-from-a-surface-for-3d-printing
% http://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-filename--varargin-

BBspectrum_norm = surf2solid(x, y, z, 'Elevation', -1);
stlwrite('BBspectrum_norm.stl',BBspectrum_norm)