###################################################################
# colorLib Maple library for use with Spectrum.mws Maple worksheet.
# Version 1.0 of 2/5/2007. Copyright (c) I.N. Galidakis, 2007.
# spectrum2xyz code based on John Walker's C code, at:
# http://www.fourmilab.ch/documents/specrend/specrend.c
# http://www.fourmilab.ch/documents/specrend/
# XYZtoCorColorTemp Maple code based on Bruce Justin Lindbloom's c
# code,at:
# http://www.brucelindbloom.com/index.html?Eqn_XYZ_to_T.html
# element spectral distributions from NIST:
# http://physics.nist.gov/cgi-bin/AtData/pt?optionslist=XXT1/
###################################################################
#Constants and variables
###################################################################
#initialize rgb color
r:=0.0:g:=0.0:b:=0.0:
###################################################################
#initialize Photoshop R, G, B color
R:=0:G:=0:B:=0:
###################################################################
#available coloring systems
NTSCsystem:=0:      #NTSC
EBUsystem:=1:       #EBU (PAL/SECAM)
SMPTEsystem:=2:     #SMPTE
HDTVsystem:=3:      #HDTV
CIEsystem:=4:       #CIE
Rec709system:=5:    #CIE REC709
###################################################################
#indexes for array of sys data
xRed:=0:yRed:=1:
xGreen:=2:yGreen:=3:
xBlue:=4:yBlue:=5:
xWhite:=6:yWhite:=7:
###################################################################
#initialize system primaries for systems above
systemSpec:=array(NTSCsystem..Rec709system,xRed..yWhite,[
    [0.67,0.33,0.21,0.71,0.14,0.08,0.3101,0.3162],
    [0.64,0.33,0.29,0.60,0.15,0.06,0.3127,0.3291],
    [0.630,0.340,0.310,0.595,0.155,0.070,0.3127,0.3291],
    [0.670,0.330,0.210,0.710,0.150,0.060,0.3127,0.3291],
    [0.7355,0.2645,0.2658,0.7243,0.1669,0.0085,0.33333333,0.33333333],
    [0.64,0.33,0.30,0.60,0.15,0.06,0.3127,0.3291]
    ]):
###################################################################
#one of the above. Change this for a different coloring system
s:=CIEsystem:
###################################################################
#initialize primaries for system
sxRed:=systemSpec[s,xRed]:syRed:=systemSpec[s,yRed]:
sxGreen:=systemSpec[s,xGreen]:syGreen:=systemSpec[s,yGreen]:
sxBlue:=systemSpec[s,xBlue]:syBlue:=systemSpec[s,yBlue]:
sxWhite:=systemSpec[s,xWhite]:syWhite:=systemSpec[s,yWhite]:
###################################################################
#gamma for LCD screen. Change for other screens
sGamma:=2.5:
###################################################################
#spectrum name (string):
spectrum:="";
###################################################################
#spectrum bounds
sStart:=100: #UV
sEnd:=2000; #IR
###################################################################
#nm bins:
kMaxBins:=sEnd-sStart;
###################################################################
#reciprocal color temp (K) for black body calculations
rt:=array(0..30,[
         0.0, 10.0e-6,20.0e-6,30.0e-6,40.0e-6,50.0e-6,
         60.0e-6,70.0e-6,80.0e-6,90.0e-6,100.0e-6,125.0e-6,
        150.0e-6,175.0e-6,200.0e-6,225.0e-6,250.0e-6,275.0e-6,
        300.0e-6,325.0e-6,350.0e-6,375.0e-6,400.0e-6,425.0e-6,
        450.0e-6,475.0e-6,500.0e-6,525.0e-6,550.0e-6,575.0e-6,
        600.0e-6
        ]):
###################################################################
#uvt coordinates array
uvt:=array(0..30,0..2,[
        [0.18006, 0.26352, -0.24341],
        [0.18066, 0.26589, -0.25479],
        [0.18133, 0.26846, -0.26876],
        [0.18208, 0.27119, -0.28539],
        [0.18293, 0.27407, -0.30470],
        [0.18388, 0.27709, -0.32675],
        [0.18494, 0.28021, -0.35156],
        [0.18611, 0.28342, -0.37915],
        [0.18740, 0.28668, -0.40955],
        [0.18880, 0.28997, -0.44278],
        [0.19032, 0.29326, -0.47888],
        [0.19462, 0.30141, -0.58204],
        [0.19962, 0.30921, -0.70471],
        [0.20525, 0.31647, -0.84901],
        [0.21142, 0.32312, -1.0182],
        [0.21807, 0.32909, -1.2168],
        [0.22511, 0.33439, -1.4512],
        [0.23247, 0.33904, -1.7298],
        [0.24010, 0.34308, -2.0637],
        [0.24792, 0.34655, -2.4681],
# Note: 0.24792 is a corrected value for the error
# found in W&S as 0.24702
        [0.25591, 0.34951, -2.9641],
        [0.26400, 0.35200, -3.5814],
        [0.27218, 0.35407, -4.3633],
        [0.28039, 0.35577, -5.3762],
        [0.28863, 0.35714, -6.7262],
        [0.29685, 0.35823, -8.5955],
        [0.30505, 0.35907, -11.324],
        [0.31320, 0.35968, -15.628],
        [0.32129, 0.36011, -23.325],
        [0.32931, 0.36038, -40.770],
        [0.33724, 0.36051, -116.45]
]):
###################################################################
#type of spectrum
emSpectrum:=TRUE;
abSpectrum:=FALSE;
###################################################################
#maximum number of points for CIE diagram
kMaxCIEPoints:=100000;
###################################################################
#points for CIE chromaticity diagram
CIEpoints:=array(1..kMaxCIEPoints,1..5);
###################################################################
#initialize coordinates to 0
X:=0.0;Y:=0.0;Z:=0.0;
x:=0.0:y:=0.0:z:=0.0:
###################################################################
#Correlated Color Temperature
CCT:=-1;
###################################################################
#Full Spectrum index
FSCI:=-1;
cspdt:=0.0;
###################################################################
#grey color
macro(grey=COLOR(RGB,0.5,0.5,0.5)):
###################################################################
#spectrum lines as a list
specList:=[];
###################################################################
#wavelength intensity of line, initialized non-proceduraly
specIntens:=array(0..kMaxBins):
for n from 0 to kMaxBins do
 specIntens[n]:=0;
od:
###################################################################
# same as above for saving
saveSpecIntens:=array(0..kMaxBins):
for n from 0 to kMaxBins do
 saveSpecIntens[n]:=0;
od:
n:='n'; #just in case we need it later...
###################################################################
#r,g,b color of spectrum bins. Calculated on demand or read from file
binColor:=array(0..kMaxBins,0..3):
###################################################################
#color matching table from:
#http://cvrl.ioo.ucl.ac.uk/
#will be read from file
cieColorMatch:=array(0..kMaxBins,0..2):
###################################################################
epsilon:=1e-4;
maxLCCT:=10000; #practically infinity
maxLLife:=24*365*10;#10 years. practically infinity
dCCT:=5785;
m:=0; #number of lamps read from file
###################################################################
#daylight's coordinates
sx:=0.325512; 
sy:=0.335131;
###################################################################
#lamp data for max 30 lamps
kLamp:=0;
kCCT:=kLamp+1;
kdbb:=kCCT+1;
kdd:=kdbb+1;
kef:=kdd+1;
klif:=kef+1;
kcost:=klif+1;
kfs:=kcost+1;
LD:=array(1..30,kLamp..kfs,[
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
	["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0],
    ["",0,0,0,0,0,0,0]
    ]):
###################################################################
# procedures and functions
###################################################################
Tr:=265;#max T radiating without thermal broadening for Na (C)
c:=299792458;#m/s
k_B:=1.380649*10^(-23);#J/K
h:=6.62607015*10^(-34);#J*s
###################################################################
#linear interpolation
LERP:=(a,b,c)->((b - a) * c + a);
###################################################################
K:=C->C+273.15;#Celsius -> Kelvin
C:=K->K-273.15;#Kelvin -> Celsius
###################################################################
#bin to wavelength
bin2wavelength:=i->sStart+i;
###################################################################
#wavelength2bin
wavelength2bin:=lambda->round(lambda)-sStart;
###################################################################
#Estimate color temperature.
XYZtoCorColorTemp:=proc()
global X,Y,Z,rt,uvt;
local us,vs,p,di,dm,i,temp;
 if X < 1.0e-20 and Y < 1.0e-20 and Z < 1.0e-20 then
  RETURN(-1); #protect against possible divide-by-zero failure
 else
  us:=4.0*X/(X+15.0*Y+3.0*Z);
  vs:=6.0*Y/(X+15.0*Y+3.0*Z);
  dm:=0.0;
  for i from 0 to 30 do
   di:=vs-uvt[i,1]-uvt[i,2]*(us-uvt[i,0]);
   if ((i > 0) and (((di < 0.0) and (dm >= 0.0))
         or ((di >= 0.0) and (dm < 0.0)))) then
    break;  #found lines bounding (us, vs) : i-1 and i
   fi;
   dm:=di;
  od;
  if i=31 then
#bad XYZ input, color temp would be less than minimum of 1666.7
# degrees, or too far towards blue
   RETURN(-1);
  fi;
  di:=di/sqrt(1.0+uvt[i,2]^2);
  dm:=dm/sqrt(1.0+uvt[i-1,2]^2);
  p:=dm/(dm-di);  # p = interpolation parameter, 0.0:i-1,1.0:i
  p:=1.0/LERP(rt[i-1],rt[i],p);
  temp:=p;
  RETURN(temp);      # success
 fi;
end:
###################################################################
#returns true if no visible lines!
specVisible:=proc()
global specIntens;
local i,vis;
 vis:=FALSE;
 for i from wavelength2bin(380) to wavelength2bin(700) do
  if specIntens[i]>0 then
   vis:=TRUE;
   break;
  fi;
 od;
 RETURN(vis);
end:
###################################################################
#returns true if no lines in spectrum!
specEmpty:=proc()
global specIntens;
local i,empty;
 empty:=TRUE;
 for i from 0 to kMaxBins do
  if specIntens[i]>0 then
   empty:=FALSE;
   break;
  fi;
 od;
 RETURN(empty);
end:
###################################################################
#read color matching table from file
#with data gotten from link above
readcieColorMatch:=proc(filename)
global cieColorMatch;
local sd,i,data;
 for i from wavelength2bin(100) to wavelength2bin(379) do
  cieColorMatch[i,0]:=0.0; #UV is invisible
  cieColorMatch[i,1]:=0.0;
  cieColorMatch[i,2]:=0.0;
 od;
 sd:=fopen(cat(dir,filename),READ,TEXT);
 for i from wavelength2bin(380) to wavelength2bin(700) do
  data:=fscanf(sd,"%d %g %g %g");
  cieColorMatch[i,0]:=data[2]; #read visible color matching values
  cieColorMatch[i,1]:=data[3];
  cieColorMatch[i,2]:=data[4];
 od;
 fclose(sd);
 for i from wavelength2bin(701) to wavelength2bin(2000) do
  cieColorMatch[i,0]:=0.0; #IR is invisible
  cieColorMatch[i,1]:=0.0;
  cieColorMatch[i,2]:=0.0;
 od;
end:
###################################################################
CIEcolorMatchGraph:=proc()
global cieColorMatch;
local i,titles,labelss,low,hiw,cieR,cieG,cieB,p1,p2,p3,p;
 titles:="CIE 1931 color matching functions x( ),y( ),z( )";
 labelss:=["nm","wgt"];
 cieR:=[];
 cieG:=[];
 cieB:=[];
 low:=380;
 hiw:=700;
 for i from wavelength2bin(low) to wavelength2bin(hiw) do
  cieR:=[op(cieR),[bin2wavelength(i),cieColorMatch[i,0]]];
  cieG:=[op(cieG),[bin2wavelength(i),cieColorMatch[i,1]]];
  cieB:=[op(cieB),[bin2wavelength(i),cieColorMatch[i,2]]];
 od;
 p1:=plot(cieR,color=red,
            axesfont=[TIMES,ROMAN,8],
            labels=labelss,
            labelfont=[TIMES,ROMAN,8],
            title=titles,
            titlefont=[TIMES,ROMAN,8]):
 p2:=plot(cieG,color=green,
            axesfont=[TIMES,ROMAN,8],
            labels=labelss,
            labelfont=[TIMES,ROMAN,8],
            title=titles,
            titlefont=[TIMES,ROMAN,8]):
 p3:=plot(cieB,color=blue,
            axesfont=[TIMES,ROMAN,8],
            labels=labelss,
            labelfont=[TIMES,ROMAN,8],
            title=titles,
            titlefont=[TIMES,ROMAN,8]):
 p:={p1} union {p2} union {p3};
 RETURN(p);
end:
###################################################################
#read all line rgb color values from file
#to save calculation time
readColors:=proc(filename)
global binColor;
local sd,i,lo,hi,data;
 lo:=wavelength2bin(100);
 hi:=wavelength2bin(379);
 for i from lo to hi do
  binColor[i,0]:=0.5; #black for UV
  binColor[i,1]:=0.5;
  binColor[i,2]:=0.5;
 od;
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 sd:=fopen(cat(dir,filename),READ,TEXT);
 for i from lo to hi do
  data:=fscanf(sd,"%d %g %g %g");
  binColor[i,0]:=data[2]; #read colors for visible
  binColor[i,1]:=data[3];
  binColor[i,2]:=data[4];
 od;
 fclose(sd);
 lo:=wavelength2bin(701);
 hi:=wavelength2bin(2000);
 for i from lo to hi do
  binColor[i,0]:=0.5; #black for IR
  binColor[i,1]:=0.5;
  binColor[i,2]:=0.5;
 od;
end:
###################################################################
#clear spectrum parameters
clearSpectrum:=proc()
global specList,specIntens,kMaxBins,
lo,hi,x,y,z,X,Y,Z,r,g,b,R,G,B,spectrum;
local i;
 if nargs=0 then
  lo:=0;
  hi:=kMaxBins;
 else
  lo:=args[1];
  hi:=args[2];
 fi;
 spectrum:="";
 x:=0.0;y:=0.0;z:=0.0;
 X:=0.0;Y:=0.0;Z:=0.0;
 r:=0.0;g:=0.0;b:=0.0;
 R:=0;G:=0;B:=0;
 specList:=[];
 for i from lo to hi do
  specIntens[i]:=0;
 od;
end:
###################################################################
#calculate CIE X,Y,Z coordinates and then
#chromaticity coordinates x,y,z
spectrum2xyz:=proc()
global x,y,z,X,Y,Z,cieColorMatch,specIntens;
local i;
 for i from wavelength2bin(380) to wavelength2bin(700) do
  X:=X+specIntens[i]*cieColorMatch[i,0];
  Y:=Y+specIntens[i]*cieColorMatch[i,1];
  Z:=Z+specIntens[i]*cieColorMatch[i,2];
 od;
 x:=X/(X+Y+Z);
 y:=Y/(X+Y+Z);
 z:=1-(x+y);
end:
###################################################################
#load line and intensity for element with participating percentage
loadLine:=proc(emission,lambda,intens,ratio)
global specIntens;
local i;
 i:=wavelength2bin(lambda);#find bin
 if emission=TRUE then
  specIntens[i]:=specIntens[i]+ratio*intens;
 else
  specIntens[i]:=specIntens[i]-ratio*intens;
 fi;
end:
###################################################################
#show distribution (used in debugging only)
showLines:=proc()
global specIntens,kMaxBins;
local i,lambda,lo,hi;
 if nargs=0 then
  lo:=0;
  hi:=kMaxBins;
 else
  lo:=wavelength2bin(args[1]);
  hi:=wavelength2bin(args[2]);
 fi;
 for i from lo to hi do
  lambda:=bin2wavelength(i);
  print(lambda,specIntens[i]);
 od;
end:
###################################################################
#x,y,z chromaticity to r,g,b conversion
xyz2rgb:=proc(xc,yc,zc)
global xRed,yRed,xGreen,yGreen,xBlue,yBlue,xWhite,yWhite,r,g,b;
local xr,yr,zr,xg,yg,zg,xb,yb,zb,xw,yw,zw,rx,ry,rz,gx,gy,gz,bx,By,
bz,rw,gw,bw;
 xr:= sxRed;yr:=syRed;zr:=1-(xr + yr);
 xg:=sxGreen;yg:=syGreen;zg:=1-(xg + yg);
 xb:=sxBlue;yb:=syBlue;zb:=1-(xb + yb);
 xw:=sxWhite;yw:=syWhite;zw:=1-(xw + yw);
 rx:=(yg*zb)-(yb*zg);ry:=(xb*zg)-(xg*zb);rz:=(xg*yb)-(xb*yg);
 gx:=(yb*zr)-(yr*zb);gy:=(xr*zb)-(xb*zr);gz:=(xb*yr)-(xr*yb);
 bx:=(yr*zg)-(yg*zr);By:=(xg*zr)-(xr*zg);bz:=(xr*yg)-(xg*yr);
 rw:=((rx*xw)+(ry*yw)+(rz*zw))/yw;
 gw:=((gx*xw)+(gy*yw)+(gz*zw))/yw;
 bw:=((bx*xw)+(By*yw)+(bz*zw))/yw;
 rx:=rx/rw;ry:=ry/rw;rz:=rz/rw;
 gx:=gx/gw;gy:=gy/gw;gz:=gz/gw;
 bx:=bx/bw;By:=By/bw;bz:=bz/bw;
 r:=(rx*xc)+(ry*yc)+(rz*zc);
 g:=(gx*xc)+(gy*yc)+(gz*zc);
 b:=(bx*xc)+(By*yc)+(bz*zc);
end:
###################################################################
#r,g,b to r,g,b inside GAMUT triangle
constrainrgb:=proc(oldr,oldg,oldb)
global r,g,b;
local w;
 if (oldr>0) then w:=0;
 else w:=oldr;fi;
 if (oldg<=w) then w:=oldg;fi;
 if (oldb<=w) then w:=oldb;fi;
 w:=-w;
 if w>0 then
  r:=oldr+w;
  g:=oldg+w;
  b:=oldb+w;
 fi;
end:
###################################################################
#gamma correction for r,g,b
gammaCorrect:=proc(c)
global sGamma;
local gamma,cc,newc;
 gamma:=sGamma;
 cc:=0.018;
 if (gamma=0) then
  if (c < cc) then
   newc:=((1.099*cc^0.45)-0.099)/cc;
  else
   newc:=(1.099*c^0.45)-0.099;
 fi;
 else
  newc:=c^(1/gamma);
 fi;
 RETURN(newc);
end:
###################################################################
#gamma correction for all 3 r,g,b values
gammaCorrectrgb:=proc(oldr,oldg,oldb)
global r,g,b;
 r:=gammaCorrect(oldr);
 g:=gammaCorrect(oldg);
 b:=gammaCorrect(oldb);
end:
###################################################################
#normalize r,g,b values
normalizergb:=proc(oldr,oldg,oldb)
global r,g,b;
local greatest;
 greatest:=max(oldr,oldg,oldb);
 if (greatest>0.0) then
  r:=oldr/greatest;
  g:=oldg/greatest;
  b:=oldb/greatest;
 fi;
end:
###################################################################
#r,g,b to Photoshop R,G,B
rgb2PhotoshopRGB:=proc(r,g,b)
global R,G,B;
 R:=round(r*255);
 G:=round(g*255);
 B:=round(b*255);
end:
###################################################################
#main processing. Calls above
calculateParameters:=proc()
global X,Y,Z,CCT;
 if specVisible()=TRUE then
  print("spectrum: ",spectrum);
  spectrum2xyz();print("chromaticity: ",x,y,z);
  CCT:=XYZtoCorColorTemp();
  if CCT=-1 then
   print("CCT not available");
  else
   print("CCT: ",round(CCT), " (K)");
  fi;
  xyz2rgb(x,y,z);#print("original rgb: ",r,g,b);
  if (r<0) or (g<0) or (b<0) then
   constrainrgb(r,g,b);#print("constrained rgb: ",r,g,b);
  fi;
  gammaCorrectrgb(r,g,b);#print("gamma corrected rgb: ",r,g,b);
  normalizergb(r,g,b);#print("normalized rgb: ",r,g,b);
  rgb2PhotoshopRGB(r,g,b);print("Photoshop RGB: ",R,G,B);
 else
  ERROR("Spectrum is empty!");
 fi;
end:
###################################################################
#spectral distribution
SPD:=proc(lambda)
global specIntens,cspdt;
local lo,hi,lambdabin,i;
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 lambdabin:=wavelength2bin(lambda);
 RETURN(specIntens[lambdabin]/cspdt);
end:
###################################################################
SPDp:=proc(lambda)
if 380<=lambda and lambda <700 then
 SPD(lambda);
elif 700<=lambda and lambda<=320+700 then
 SPD(lambda-320);
else
 -10000;#for debuging
fi;
end:
###################################################################
#cumulative Equal Energy distribution
CEE:=proc(lambda)
 RETURN(evalf((lambda-380)/(700-380)));
end:
###################################################################
#calculate FSCI
FSi:=proc()
global cspdt,FSCI,specIntens;
local lo,hi,lambda,lambdac,w,dw,cSPD,s,i,fs,fsi,fsci;
 s:=0.0;dw:=10;
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 cspdt:=add(specIntens[i],i=lo..hi);
 for w from 0 to 700-380 by dw do
  for lambdac from 380+w to 700+w do
   cSPD:=add(SPDp(lambda),lambda=380+w..lambdac);
   s:=s+(cSPD-CEE(lambdac-w))^2;
  od;
 od;
 fsi:=s*dw/(700-380);
 fsci:=100-5.1*fsi;
 if fsci<=0 then
  FSCI:=0;
 else
  FSCI:=fsci;
 fi;#print("Full spectrum index:", FSCI);
end:
###################################################################
#fade colors in 380-399. Linear fade, simulates darker purples
fadeUV:=proc()
global binColor,cieColorMatch;
local i,lo,hi,mlo,mhi;
 for i from wavelength2bin(380) to wavelength2bin(380) do
  binColor[i,0]:=(i-wavelength2bin(380))*binColor[i,0]
    /(wavelength2bin(400)-wavelength2bin(380));
  binColor[i,1]:=(i-wavelength2bin(380))*binColor[i,1]
    /(wavelength2bin(400)-wavelength2bin(380));
  binColor[i,2]:=(i-wavelength2bin(380))*binColor[i,2]
    /(wavelength2bin(400)-wavelength2bin(380));
 od;
end:
###################################################################
#fade colors in 675-700. Linear fade, simulates darker reds
fadeIR:=proc()
global binColor,cieColorMatch;
local i;
 for i from wavelength2bin(651) to wavelength2bin(700) do
  binColor[i,0]:=(wavelength2bin(700)-i)*binColor[i,0]
    /(wavelength2bin(700)-wavelength2bin(650));
  binColor[i,1]:=(wavelength2bin(700)-i)*binColor[i,1]
    /(wavelength2bin(700)-wavelength2bin(650));
  binColor[i,2]:=(wavelength2bin(700)-i)*binColor[i,2]
    /(wavelength2bin(700)-wavelength2bin(650));
 od;
end:
###################################################################
#Calculate pure colors for wavelengths. Values have been stored in
#file 'colors.txt', so not used unless a recalculation is called for
calculateBinColors:=proc()
global binColor,specIntens,x,y,z,r,g,b;
local i,lo,hi;
 lo:=wavelength2bin(100);
 hi:=wavelength2bin(379);
 for i from lo to hi do
  binColor[i,0]:=0.5; #UV is invisible
  binColor[i,1]:=0.5;
  binColor[i,2]:=0.5;
 od;
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 for i from lo to hi do
  clearSpectrum(lo,hi);
  specIntens[i]:=1000;
  spectrum2xyz();
  xyz2rgb(x,y,z);
  if (r<0) or (g<0) or (b<0) then
   constrainrgb(r,g,b);
  fi;
  gammaCorrectrgb(r,g,b);
  normalizergb(r,g,b);
  binColor[i,0]:=r;
  binColor[i,1]:=g;
  binColor[i,2]:=b;
 od;
 lo:=wavelength2bin(701);
 hi:=wavelength2bin(2000);
 for i from lo to hi do
  binColor[i,0]:=0.5; #IR is invisible
  binColor[i,1]:=0.5;
  binColor[i,2]:=0.5;
 od;
 fadeUV();
 fadeIR();
end:
###################################################################
#normalize line intensities. Max=100
normalizeIntens:=proc()
global specIntens;
local i,maxIntens;
 maxIntens:=0;
 for i from 0 to kMaxBins do
  maxIntens:=max(maxIntens,specIntens[i]);
 od;
 for i from 0 to kMaxBins do
  specIntens[i]:=specIntens[i]/maxIntens*100;
 od;
end:
###################################################################
#create lines list
createList:=proc()
global specIntens,specList;
local i;
 for i from 0 to kMaxBins do
  if specIntens[i]>0 then
   specList:=[op(specList),[i,specIntens[i]]];
  fi;
 od;
end:
###################################################################
#render spectrum.
#lo=lower bound in nm >=100
#hi=upper bound in nm <=2000
#maxIntensShow=intensity as a percentage, eg, 50,100, etc
renderSpectrum:=proc(lo,hi,maxIntensShow,inverted)
global binColor,specList,spectrum,emptySpec;
local i,j,p,los,his,mr,mg,mb,ss;
 if specEmpty()=FALSE then
  normalizeIntens();
  createList();
  ss:={};
  los:=max(lo,100);
  his:=min(2000,hi);
  for j from 1 to nops(specList) do
   i:=specList[j,1];
   if (los<=bin2wavelength(i)) and (bin2wavelength(i)<=his) then
    if inverted=TRUE then
     if 380<=bin2wavelength(i) and bin2wavelength(i)<=700 then
      mr:=1-binColor[i,0];
      mg:=1-binColor[i,1];
      mb:=1-binColor[i,2];
     else
      mr:=binColor[i,0]; #UV/IR
      mg:=binColor[i,1];
      mb:=binColor[i,2];
     fi;
    else #not inverted
     mr:=binColor[i,0];
     mg:=binColor[i,1];
     mb:=binColor[i,2];
    fi;
    p:=plot([[bin2wavelength(i),0],[bin2wavelength(i),specIntens[i]]],
      los..his,
      0..maxIntensShow,
      axesfont=[TIMES,ROMAN,8],
      labels=["nm","R.I."],
      labelfont=[TIMES,ROMAN,8],
      color=COLOR(RGB,mr,mg,mb),
      title=spectrum,
      titlefont=[TIMES,ROMAN,9]
        );
    ss:=ss union {p};
   fi;
  od;
  display(ss,size=[393,404]);
 else
  ERROR("Spectrum is empty!");
 fi;
end:
###################################################################
#loads specific element from file, with strength ratio
###################################################################
#loads specific element from file, with strength ratio
loadElement:=proc(emission,element,ratio)
global spectrum;
local n,sd,data,lambda,intens,pm;
 sd:=fopen(cat(dir,element, ".txt"),READ,TEXT);
 print("element: ",element);
 #print("ratio: ",ratio);
 for n from 1 to 100000 do
  data:=fscanf(sd,"%g %g");
  if feof(sd) then break
  else
   lambda:=data[1];
   intens:=data[2];
   loadLine(emission,lambda,intens,ratio);
  fi;
 od:
 fclose(sd);
 if emission=TRUE then
  pm:="+";
 else
  pm:="-";
 fi;
 if spectrum="" then
  if ratio<>1 then
   spectrum:=cat(element,"(",convert(ratio,string),")");
  else
   spectrum:=element;
  fi;
 else
  if ratio<>1 then
   spectrum:=cat(spectrum,pm,element,"(",convert(ratio,string),")");
  else
   spectrum:=cat(spectrum,pm,element);
  fi;
 fi;
 #print("lines: ",n-1);
end:
###################################################################
#load element combinations
loadElements:=proc()#arguments: emission type,X1,r1,...,Xn,rn
local i,emission,element,ratio;
 if nargs-1 mod 2<>0 then
  ERROR("Elements/ratios should come in pairs!");
 fi;
 emission:=args[1];
 for i from 2 to nargs by 2 do
  element:=args[i];
  ratio:=args[i+1];
  loadElement(emission,element,ratio);
 od;
end:
###################################################################
#append data for sources to file
appendData:=proc(newspec,eff,life,cost,filename)
global x,y,R,G,B,CCT,spectrum;
local sd;
 sd:=fopen(cat(dir,filename),APPEND,TEXT);
 fprintf(sd,"%s %g %g %g %g %g %d %d %d %d %g\n",newspec,x,y,R,G,B,round(CCT),round(eff),round(life),round(cost),FSCI);
 fclose(sd);
end:
###################################################################
#calculate intensity of black body of temperature=T at lambda
blackBody:=proc(T,lambda)
local lambdam;
 lambdam:=lambda/10^9;#nm 2 m
 8*Pi*h*c/lambdam^5/(exp(h*c/(lambdam*k_B*T))-1);
end:
###################################################################
#load black body of temp T for wavelengths low-hiw (nm)
loadBlackBody:=proc(T,low,hiw)
global spectrum,specIntens;
local i,maxIntens,lambda,intens,lo,hi;
 spectrum:=cat(spectrum,"+Black Body @",convert(T,string),"K");
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 maxIntens:=0.0;
 for i from lo to hi do #find maximum for range
  lambda:=bin2wavelength(i);
  maxIntens:=max(maxIntens,blackBody(T,lambda));
 od;
 for i from lo to hi do
  lambda:=bin2wavelength(i);
  intens:=specIntens[i]+blackBody(T,lambda)/maxIntens*100;
  loadLine(emSpectrum,lambda,intens,1);
 od;
end:
###################################################################
#write all rgb values to file
#used to store rgb values after first calculation
#to have them handy later
writeColors:=proc(filename)
global binColor;
local sd,i,lo,hi;
 sd:=fopen(cat(dir,filename),APPEND,TEXT);
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 for i from lo to hi do
  fprintf(sd,"%d %10.8g %10.8g %10.8g\n",bin2wavelength(i),
        binColor[i,0],binColor[i,1],binColor[i,2]);
 od;
 fclose(sd);
end:
###################################################################
#transformations from CIE 1960 to CIE 1931
xT:=(u,v)->3*u/(2*u-8*v+4);
yT:=(u,v)->2*v/(2*u-8*v+4);
###################################################################
#slope at (x,y) on plankian locus
tT:=(x,y,t)->(2*y*(4*t-1)+2*t)/(2*x*(4*t-1)+3);
###################################################################
#a general line
linef:=(m,x0,y0,x)->m*(x-x0)+y0;
###################################################################
#calculate Judd isothermal
juddLine:=proc(i)
global uvt;
local txt,T,m,x0,y0,x,dx,tp,aligntp,p1,p2,p;
x0:=xT(uvt[i,0],uvt[i,1]);
y0:=yT(uvt[i,0],uvt[i,1]);
m:=tT(x0,y0,uvt[i,2]);
dx:=cos(arctan(m))/15;
p1:=plot(linef(m,x0,y0,x),x=x0-dx..x0+dx,
        axesfont=[TIMES,ROMAN,8],
        labels=["cx","cy"],
        labelfont=[TIMES,ROMAN,8],
        color=grey,
        scaling=CONSTRAINED);
if rt[i]=0 then
 txt:="oo";
else
 T:=1/rt[i];
 if T>=10000 or round(round(T/100)/10)=round(T/100)/10 then
  txt:=convert(round(T/1000),string);
 else
  txt:=convert(evalf(round(T/100)/10,2),string);
 fi;
 
fi;
if m>1 then
 tp:=[x0+dx,linef(m,x0,y0,x0+dx),txt];
 aligntp:={ABOVE,RIGHT}
else
 tp:=[x0-dx,linef(m,x0,y0,x0-dx),txt];
 aligntp:={ABOVE,LEFT}
fi;
p2:=textplot(tp,
         axesfont=[TIMES,ROMAN,8],
         labels=["cx","cy"],
         labelfont=[TIMES,ROMAN,8],
         color=grey,
         font=[HELVETICA,6],
         align=aligntp,
         scaling=CONSTRAINED);
p:={p1} union {p2};
RETURN(p);
end:
###################################################################
#generate a graph of the planckian locus
planckianLocus:=proc()
global xyt;
local PL,p1,p2,p,i;
 print("Tabulating Planckian locus and isothermal lines...");
 PL:=[];
 p1:={};
 for i from 0 to 30 do #calculate Plankian locus and judd lines
  PL:=[op(PL),[xT(uvt[i,0],uvt[i,1]),yT(uvt[i,0],uvt[i,1])]];
  p1:=p1 union juddLine(i);
 od:
 p2:=pointplot(PL,view=[0..0.8,0..0.9],
          axesfont=[TIMES,ROMAN,8],
          labels=["cx","cy"],
          labelfont=[TIMES,ROMAN,8],
          color=grey,
          style=LINE,
          thickness=2,
          scaling=CONSTRAINED);
 p:=p1 union {p2};
 RETURN(p);
end:
###################################################################
#generate spectrum CIE chromaticity point for spectrum
lampChromaticities:=proc(filename)
local lx,ly,sd,i,lamp,data,p1,p2,p;
 p:={};
 sd:=fopen(cat(dir,filename),READ,TEXT); #load all lamps
 for i from 1 to 10000 do
  data:=fscanf(sd,"%s %g %g %d %d %d %d %d %d %d %g/n");
  if feof(sd) then break
  else
   lamp:=data[1];
   lx:=data[2];
   ly:=data[3];
   p1:=plot([[lx,ly]],0..0.8,0..0.9,
         axesfont=[TIMES,ROMAN,8],
         style=POINT,
         labels=["cx","cy"],
         labelfont=[TIMES,ROMAN,8],
         symbol=CROSS,
         color=grey,
         scaling=CONSTRAINED);
   p2:=textplot([lx,ly,lamp],
         axesfont=[TIMES,ROMAN,8],
         labels=["cx","cy"],
         labelfont=[TIMES,ROMAN,8],
         color=grey,
         align={ABOVE,RIGHT},
         scaling=CONSTRAINED);
   p:=p union {p1} union {p2};
  fi;
 od:
 fclose(sd);
 RETURN(p);
end:
###################################################################
#compromization optimization
F:=(CCT,dbb,dd,ef,lif,cost,fs,w1,w2,w3,w4,w5,w6,w7)
->w1*CCT+w2*dbb+w3*dd+w4*ef+w5*lif+w6*cost+w7*fs;
###################################################################
#normalize data
normalizeData:=proc(m)
global LD;
local i,maxCCT,maxdbb,maxdd,maxef,maxlif,maxcost,maxfsi;
 maxCCT:=-1e5;
 maxdbb:=-1e5;
 maxdd:=-1e5;
 maxef:=-1e5;
 maxlif:=-1e5;
 maxcost:=-1e5;
 maxfsi:=-1e5;
 for i from 1 to m do
  maxCCT:=max(maxCCT,LD[i,kCCT]);
  maxdbb:=max(maxdbb,LD[i,kdbb]);
  maxdd:=max(maxdd,LD[i,kdd]);
  maxef:=max(maxef,LD[i,kef]);
  maxlif:=max(maxlif,LD[i,klif]);
  maxcost:=max(maxcost,LD[i,kcost]);
  maxfsi:=max(maxfsi,LD[i,kfs]);
 od;
 for i from 1 to m do
  print("---------------------------");
  print(LD[i,kLamp]);
  LD[i,kCCT]:=100.*(1-abs((LD[i,kCCT]-dCCT)/(maxCCT-dCCT)));#print("PCCT",LD[i,kCCT]);
  LD[i,kdbb]:=100.*(1-LD[i,kdbb]/maxdbb);#print("PLB",LD[i,kdbb]);
  LD[i,kdd]:=100.*(1-LD[i,kdd]/maxdd);#print("PLD",LD[i,kdd]);
  LD[i,kef]:=100.*LD[i,kef]/maxef;#print("PLE",LD[i,kef]);
  LD[i,klif]:=100.*LD[i,klif]/maxlif;#print("PLT",LD[i,klif]);
  LD[i,kcost]:=100.*(1-LD[i,kcost]/maxcost);#print("PLC",LD[i,kcost]);
  LD[i,kfs]:=100*LD[i,kfs]/maxfsi;#print("PFS:",LD[i,kfs]);
 od;
end:
###################################################################
#sorting function for sort
S:=proc(a,b)
RETURN(evalb(a>=b));
end:
###################################################################
#find rating
rating:=proc(v,L)
local i,r;
 for i from 1 to nops(L) do
  if v=L[i] then
   r:=i;
   break;
  fi;
 od;
 RETURN(r);
end:
###################################################################
#calculate optimal lamp
optimalData:=proc(w1,w2,w3,w4,w5,w6,w7)
global LD,m;
local i,oe,oes,soes,maxoe,lp,rp,r;
 oes:=[];
 soes:=[];
 maxoe:=0;
 lp:="(";rp:=")";
 for i from 1 to m do #find maximum
  maxoe:=max(maxoe,F(LD[i,kCCT],LD[i,kdbb],LD[i,kdd],LD[i,kef],LD[i,klif],LD[i,kcost],LD[i,kfs],w1,w2,w3,w4,w5,w6,w7));
 od;
 for i from 1 to m do #calculate winners
  oe:=F(LD[i,kCCT],LD[i,kdbb],LD[i,kdd],LD[i,kef],LD[i,klif],LD[i,kcost],LD[i,kfs],w1,w2,w3,w4,w5,w6,w7);
  r:=evalf(100*oe/maxoe,3);
  oes:=[op(oes),r];
 od;
 soes:=sort(oes,S);
 for i from 1 to m do
  printf("%6s %f %s%d%s\n",LD[i,kLamp],oes[i],lp,rating(oes[i],soes),rp);
 od;
end:
###################################################################
printLampData:=proc()
global LD,m;
local i;
 printf("%6s %5s %5s %5s %5s %5s %5s %5s\n","Lamp","CCT","LB","LD","LE","LT","LC","FS");
 for i from 1 to m do
  printf("%6s %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f\n",
  LD[i,kLamp],LD[i,kCCT],LD[i,kdbb],LD[i,kdd],LD[i,kef],LD[i,klif],LD[i,kcost],LD[i,kfs]);
 od;
end:
###################################################################
#generate comparison data
calculateLampData:=proc(filename)
global m,LD,x,y,z,specIntens,lampData;
local lo,hi,low,hiw,lx,ly,jx,jy,bbx,bby,lT,sd,i,lef,llif,lcost,lfs,
data,lamp,dlbb,dld;
 clearLampData();
 low:=380;
 hiw:=700;
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 print("Calculating lamp data...");
 sd:=fopen(cat(dir,filename),READ,TEXT); #load all lamps
 for i from 1 to 10000 do
  data:=fscanf(sd,"%s %g %g %d %d %d %d %d %d %d %g/n");
  if feof(sd) then break
  else
   lamp:=data[1];
   lx:=data[2];
   ly:=data[3];
   lT:=data[7];
   lef:=data[8];
   llif:=data[9];
   lcost:=data[10];
   lfs:=data[11];
   
   LD[i,kLamp]:=lamp;
   LD[i,kCCT]:=lT;
   LD[i,kef]:=lef;
   LD[i,klif]:=llif;
   LD[i,kcost]:=lcost;
   LD[i,kfs]:=lfs;
	
	 if lT=-1 then #if no CCT
    LD[i,kCCT]:=maxLCCT; #set to large
   fi;
   
   if llif=-1 then #daylight has infinite lifetime
    LD[i,klif]:=maxLLife; #set to large
   fi;
   
   clearSpectrum(lo,hi);
   loadBlackBody(LD[i,kCCT],low,hiw);
   calculateParameters();
#x,y now hold the corresponding black body's coordinates
   bbx:=x;
   bby:=y;
   
   dlbb:=sqrt((bbx-lx)^2+(bby-ly)^2); #distance from black body
   dld:=sqrt((sx-lx)^2+(sy-ly)^2); #distance from daylight

   LD[i,kdbb]:=dlbb;
   LD[i,kdd]:=dld;
   print(cat(lamp,"..."));
  fi;
 od:
 fclose(sd);
 m:=i-1; #number of lamps;
 normalizeData(m);#normalize all data
end:
###################################################################
#generate spectrum CIE chromaticity point for spectrum
spectrumChromaticity:=proc()
global x,y,z,spectrum;
local p1,p2,p,low,hiw,lo,hi,maxIntensShow,emission;
 low:=380;
 hiw:=700;
 maxIntensShow:=100;
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 emission:=args[1];
 clearSpectrum(lo,hi);
 loadElements(emission,args[2..nargs]);
 calculateParameters();
 p:={};
 p1:=plot([[x,y]],0..0.8,0..0.9,
         axesfont=[TIMES,ROMAN,8],
         style=POINT,
         labels=["cx","cy"],
         labelfont=[TIMES,ROMAN,8],
         symbol=CROSS,
         color=grey,
         scaling=CONSTRAINED);
 p2:=textplot([x,y,spectrum],
         axesfont=[TIMES,ROMAN,8],
         labels=["cx","cy"],
         labelfont=[TIMES,ROMAN,8],
         color=grey,
         align={ABOVE,RIGHT},
         scaling=CONSTRAINED);
 p:={p1} union {p2};
 RETURN(p);
end:
###################################################################
#generate element points on CIE diagram
elementChromaticities:=proc(filename)
global x,y;
local i,sd,element,data,lo,hi,p;
 print("Tabulating element data...");
 lo:=wavelength2bin(380);
 hi:=wavelength2bin(700);
 p:={};
 sd:=fopen(cat(dir,filename),READ,TEXT); #load all elements
 for i from 1 to 10000 do
  data:=fscanf(sd,"%s");
  if feof(sd) then break
  else
   element:=data[1];
   p:=p union spectrumChromaticity(emSpectrum,element,1);
  fi;
 od:
 fclose(sd);
 RETURN(p);
end:
###################################################################
#TRUE if point x,y lies inside the poly (list). Boundary included.
pointInPolygon:=proc(N,poly,x,y)
local oddNodes,i,j;
 oddNodes:=0;
 j:=0;
 for i from 0 to N-1 do
  j:=j+1 mod N;
  if poly[i+1,2]<y and poly[j+1,2]>=y
        or poly[j+1,2]<y and poly[i+1,2]>=y then
   if poly[i+1,1]+(y-poly[i+1,2])/(poly[j+1,2]-poly[i+1,2])
          *(poly[j+1,1]-poly[i+1,1])<x then
    oddNodes:=oddNodes+1 mod 2;
   fi;
  fi;
 od;
 RETURN(oddNodes);
end:
###################################################################
#initialize CIE diagram points to save memory
initializeCIEpoints:=proc(numPoints)
global CIEpoints;
local i;
 for i from 1 to numPoints do
  CIEpoints[i,1]:=0.0;
  CIEpoints[i,2]:=0.0;
  CIEpoints[i,3]:=0.0;
  CIEpoints[i,4]:=0.0;
  CIEpoints[i,5]:=0.0;
 od;
end:

###################################################################
#create chromaticity diagram, with resolution dx
#max recommended dx=0.006. dx>0.006 takes > 20 minutes
CIEdiagram:=proc(dx,inverted)
global x,y,z,r,g,b,specIntens,CIEpoints,kMaxCIEPoints;
local bound,i,di,N,np,low,hiw,lo,hi,xx,maxx,maxy,yy,ub,cd,dr,p,CIEPL;
 print("Initializing CIE points...");
 np:=floor(0.8*0.9/dx^2); #calculate total number of points
 if np>kMaxCIEPoints then
  ERROR("Too many points in CIE diagram. This would take forever!");
  RETURN;
 fi;
 initializeCIEpoints(np);
 bound:=[];
 print("Calculating boundary...");
 low:=400;#start and end of wavelengths for CIE diagram
 hiw:=650;
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 di:=wavelength2bin(low+5)-wavelength2bin(low); #boundary resolution
 N:=0;
 #calculate bounding polygon with resolution dx
 maxx:=0.0;maxy:=0.0;
 for i from lo to hi by di do
  clearSpectrum(lo,hi);
  specIntens[i]:=1000;
  spectrum2xyz();
  maxx:=max(maxx,x);#find maximum values
  maxy:=max(maxy,y);
  N:=N+1; #new poly point
  bound:=[op(bound),[x,y]]; #build boundary
 od;
 print("Calculating list of points...");
 np:=0;
 for yy from maxy to 0 by -dx do #draw colors
  for xx from 0 to maxx by dx do
   if pointInPolygon(N,bound,xx,yy)=1 then #precalculate x,y,z,r,g,b
    x:=xx;
    y:=yy;
    z:=1-(xx+yy);
    xyz2rgb(x,y,z);
    if (r<0) or (g<0) or (b<0) then
     constrainrgb(r,g,b);
    fi;
    gammaCorrectrgb(r,g,b);
    normalizergb(r,g,b);
    np:=np+1;
    if np>kMaxCIEPoints then
     print(kMaxCIEPoints, " points reached");
     break;
    fi;
    CIEpoints[np,1]:=x;
    CIEpoints[np,2]:=y;
    if inverted=TRUE then
     CIEpoints[np,3]:=1-r;
     CIEpoints[np,4]:=1-g;
     CIEpoints[np,5]:=1-b;
    else
     CIEpoints[np,3]:=r;
     CIEpoints[np,4]:=g;
     CIEpoints[np,5]:=b;
    fi;
   fi;
  od;
  if np>kMaxCIEPoints then
     break;
  fi;
 od;
 #build graph of CIE diagram
 print("Building graph...");
 #cd:={polygon(bound)}; if we want the bounding polygon
 cd:={};
 dr:=dx/4;
 ub:=min(kMaxCIEPoints,np);
 CIEPL:=[seq([CIEpoints[i,1],CIEpoints[i,2],
                C(CIEpoints[i,3],CIEpoints[i,4],CIEpoints[i,5])],i=1..ub)];
 for i from 1 to ub do
  p:=pointplot([CIEpoints[i,1],CIEpoints[i,2]],
             axesfont=[TIMES,ROMAN,8],
             view=[0..0.8,0..0.9],
             labels=["cx","cy"],
             labelfont=[TIMES,ROMAN,8],
             color=COLOR(RGB,CIEpoints[i,3],CIEpoints[i,4],CIEpoints[i,5]),
             style=POINT,
             symbol=POINT,
             scaling=CONSTRAINED);
  cd:=cd union {p}; #add point
 od;
 RETURN(cd);
end:
###################################################################
#weighting function. Doesn't work. Don't use!
w:=(x,m)->(x-550)*m+.2698e-1;
###################################################################
averageSlope:=proc(filename)
local i,sd,data,wav,lastwav,intens,lastintens,m;
 m:=0;lastwav:=0;lastintens:=0;
 sd:=fopen(cat(dir,filename,".txt"),READ,TEXT);
 for i from 1 to 10000 do
  data:=fscanf(sd,"%g %g");
  if feof(sd) then break;
  else
   wav:=data[1];
   intens:=data[2];
   m:=(m+(intens-lastintens)/(wav-lastwav))/2;
   lastwav:=wav;
   lastintens:=intens;
  fi;
 od;
 fclose(sd);
 RETURN(m);
end:
###################################################################
#average spectrum to balance distribution
averageSpectrum:=proc(filename1,filename2)
local i,data,sd1,sd2,wavelength,intens,m;
 m:=averageSlope(filename1);
 sd1:=fopen(cat(dir,filename1,".txt"),READ,TEXT);
 sd2:=fopen(cat(dir,filename2,".txt"),WRITE,TEXT);
 for i from 1 to 10000 do
  data:=fscanf(sd1,"%g %g");
  if feof(sd1) then break;
  else
   wavelength:=data[1];
   intens:=data[2];
   fprintf(sd2,"%g %g\n", wavelength, intens*.2698e-1/w(wavelength,1.1989*m));
  fi;
 od;
 fclose(sd1);
 fclose(sd2);
end:
###################################################################
#generate emission spectrum with given parameters
emissionSpectrum:=proc()#args: low,hi,X1,p1,...,Xn,pn,maxIntensShow,inverted
 local low,hiw,lo,hi,maxIntensShow,inverted;
 low:=args[1];
 hiw:=args[2];
 inverted:=args[nargs];
 maxIntensShow:=args[nargs-1];
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 clearSpectrum(lo,hi);
 loadElements(emSpectrum,args[3..nargs-2]);
 calculateParameters();
 #FSi();
 renderSpectrum(low,hiw,maxIntensShow,inverted);
end:
###################################################################
#generate Black Body spectrum with given parameters
blackBodySpectrum:=proc(T,low,hiw,maxIntensShow,inverted)
local lo,hi;
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 clearSpectrum(lo,hi);
 loadBlackBody(T,low,hiw);
 calculateParameters();
 #FSi();
 renderSpectrum(low,hiw,maxIntensShow,inverted);
end:
###################################################################
absorptionSpectrum:=proc() #args: T,low,hiw,X1,p1,...Xn,pn,maxIntensShow,inverted
local T,lo,hi,low,hiw,maxIntensShow,inverted;
 T:=args[1];
 low:=args[2];
 hiw:=args[3];
 inverted:=args[nargs];
 maxIntensShow:=args[nargs-1];
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 clearSpectrum(lo,hi);
 loadBlackBody(T,low,hiw);
 loadElements(abSpectrum,args[4..nargs-2]);
 calculateParameters();
 #FSi();
 renderSpectrum(low,hiw,maxIntensShow,inverted);
end:
###################################################################
mixedSpectrum:=proc() #args: T,low,hiw,X1,p1,...Xn,pn,maxIntensShow,inverted
local T,lo,hi,low,hiw,maxIntensShow,inverted;
 T:=args[1];
 low:=args[2];
 hiw:=args[3];
 inverted:=args[nargs];
 maxIntensShow:=args[nargs-1];
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 clearSpectrum(lo,hi);
 loadBlackBody(T,low,hiw);
 loadElements(emSpectrum,args[4..nargs-2]);
 calculateParameters();
 #FSi();
 renderSpectrum(low,hiw,maxIntensShow,inverted);
end:
###################################################################
maxIntens:=proc(cSpecIntens,lo,hi)
local maxIntens,i;
maxIntens:=0;
for i from lo to hi do
 maxIntens:=max(cSpecIntens[i],maxIntens);
od;
maxIntens;
end:
###################################################################
#TRUE if isolated emission line
isoLine:=proc(lo,hi,i)
global specIntens;
local insp,isol,strong;
insp:=(lo<i) and (i<hi);
if insp then
 isol:=(specIntens[i]>specIntens[i-1]) and (specIntens[i]>specIntens[i+1]);
else
 isol:=FALSE;
fi;
RETURN(insp and isol);
end:
###################################################################
#dirac delta for single line in low temp spectrum
delt:=(x,s)->piecewise(x=s,1,0):
###################################################################
#Lorentzian profile distribution for high T/pressure lines
Lorentzian:=(x,cLambda,gam,h,ir)->h/(1+((x-cLambda)/(3.3*ir*gam))^2);
###################################################################
LorentzianLine:=proc(x,cLambda,T1,T2,h,ir)
if T1>Tr then
 if T2>=Tr then#there is emission AND absorption absorption
  Lorentzian(x,cLambda,T1/Tr-1,h,ir)-Lorentzian(x,cLambda,T2/Tr-1,h,ir);
 else#T2<Tr, so no absorption
  Lorentzian(x,cLambda,T1/Tr-1,h,ir);
 fi;
else#T1<=Tr
 h*delt(x,cLambda);
fi;
end:
###################################################################
#load Lorentzian distribution for wavelengths low to hiw
loadLorentzians:=proc(low,hiw,T1,T2)
global spectrum,specIntens;
local i,j,lambda,intens,lo,hi,maxi,centerLambda,ir;
local cSpecIntens:=array(0..kMaxBins);
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 #print(589,specIntens[wavelength2bin(589)]);
 for i from lo to hi do #temporary load copy of intensities
  cSpecIntens[i]:=specIntens[i];
 od;
 maxi:=maxIntens(cSpecIntens,lo,hi);#precalculate maximum intensity line
 #print(maxi);
 for i from lo to hi do #clear all intensities
  specIntens[i]:=0;
 od;
 for i from lo to hi do
  lambda:=bin2wavelength(i);
	intens:=0;#start with null bin
	for j from lo to hi do
	 if cSpecIntens[j]>0 then#if original spectrum has a line there
	  centerLambda:=bin2wavelength(j);
		ir:=evalf(cSpecIntens[j]/maxi);#line intens ratio relative to strongest
	  intens:=intens+round(evalf(LorentzianLine(lambda,centerLambda,
		T1,ir*T2,cSpecIntens[j],ir)));#add Lorentzian contrib from line j
	 fi;
	od;
	specIntens[i]:=intens;
 od;
end:
###################################################################
hpSpectrum:=proc()
#(low,hiw,elf1,elp1,...elfn,elpn,T1,T2,T3,maxIntensShow,inverted)
local lo,hi,low,hiw,T1,T2,T3,maxIntensShow,inverted;
 low:=args[1];
 hiw:=args[2];
 lo:=wavelength2bin(low);
 hi:=wavelength2bin(hiw);
 inverted:=args[nargs];
 maxIntensShow:=args[nargs-1];
 T3:=args[nargs-2];
 T2:=args[nargs-3];
 T1:=args[nargs-4];
 clearSpectrum(lo,hi);
 loadElements(emSpectrum,args[3..nargs-5]);
 loadLorentzians(low,hiw,T1,T2);
 loadBlackBody(T3,low,hiw);
 #showLines(500,530);
 calculateParameters();
 #FSi();
 renderSpectrum(low,hiw,maxIntensShow,inverted);
end:
###################################################################
# Code starts here
###################################################################
#initialize colors and color match
readcieColorMatch("CIE1931colorMatch.txt"):
###################################################################
#to calculate all bin colors from scratch, erase file "colors.txt"
#and uncomment the next two lines to write the new colors. Then
#comment them back
#calculateBinColors():
#writeColors("colors.txt"):
###################################################################
readColors("colors.txt"):
