#Dispersion.rb implements the gaussian dispersion model using the ruby language
#Bugs,Comments or Questions unomi@techemail.com


require "opengl"
require "glut"
require "mathn"


#Hi and welcome to dispersion
#
#This program is written in a language called ruby
#im going to have to show you a bit of this gem 
#for you to ba able to alter and appreciate this program

#I am pretty certain that reading thru this program and acquainting
#yourself with the ruby language will be delicious

#this program will open files, readlines,split translate each in turn

#do not become alarmed =) ALSO if you are about to stop reading for any reason, hold it off 20 secs, 
#and start reading upwards from doModel() near the bottom 

#programs basically do stuff to stuff,
#stuff they do stuff to they call variables

#ruby is pretty nice with variables
E = 1 # emission rate g m-1s-1

#now we have defined a variable called E that holds a value of 1

#ruby has alot of types of variables or only like 1
#either way 

#variable names are pretty open, you dont have to define something as a number or a string or whatever if you dont want

#the $ sign means that its a global variable so we can get to it from inside functions/methods (google define:scope)
#with E above, the fact that the first letter is in uppercase means that it is a Constant, it cannot be changed during the run
#of the program

#sometimes you want to store complex information and you use more refined classes of variables
$poi = Array.new #defines a point of interest as a list/array
$poi[0]= 200 #x value
 $poi[2] = 4 # z value

#an array is basically a list, starting at 0, so the first value of the variable $poi 
#is 200  and the 3rd is 4, reason being anticipation of doing it fully 3d

$abcdf = Hash.new 
#A Hash is like a list but instead of numbers they have names like 'larsen' or 'F' 


#variables also like pointing to files
 $metf= File.new("met.txt")  #associates $metf with met.txt
#the met.txt has the metereological data in time slices of an hour  

#this jumps 200 timeslices forward in the file
#loops are quite nice and this is one of them
#the |h| thing isnt used, but would hold 0,1,2,3,4.. and so on to 200 
# for example putting 
#print h, $metf.readline 
#in the loop would print out the first 200 
#lines prefixed with the line number
200.times do |h| 
	  $metf.readline
  end


#programs also need functions to be really cool
#you can have programs with no functions, but they arent cool, ask anyone.
#in ruby you define a function like youll see in the first example =)

#the first example will be sigmaZbalance
#its a function that reads in information from a file 
#if you dont know what a loop is or any other questions, contact me asap. 
#then just get a 'feel' for what is going on



def sigmaZbalance() #define function sigmaZbalance
# reads in the weights for computing sigmaZ based on Pasquill atmospheric stability guide Ep1 table 1  
  abcdfile = File.new("abcdf.txt") #open it and see for yourself 

	('A'..'F').each do |stb|  #from the letters A to F, for each of them as stb
	  $abcdf[stb]= abcdfile.readline.split   #reads in a line and splits it up according to whitespace, then puts in that list of numbers under 'A' in variable $abcdf
   end  
 
#from the letters A to F, for each of them get their list of values and convert them to numbers
  ('A'..'F').each do |stb|  #from the letters A to F, for each of them as stb
      hold=  $abcdf[stb].collect {|i|  i.to_f} #get $abcdf's list of values and convert them to numbers
      $abcdf[stb]= hold #fudge
  end 

#thats it for sigmaZbalance()
#all it does is prime or intialize $abdcf with the values from the text file.
#these values are then looked up later from 
#  getSigmaZ(x)
end

#so? how did that feel?
#what did we end up doing?
#can you think of why we did it? =)

#ok well dont worry about it too much if you have no idea what happened, just read on and come back to it later
#this program is ours, add comments where you didnt understand what was going on before
#start

$step=3
#controls how fine a resolution to use 
#when calculating concentration maps
#can make everything very slow 1=slowest 5=default

$cutoff=0.0001 #completely spurious numbers used for deciding what concentration ranges are interesting		
 $climit=0.0002
 $c=0 
 

sigmaZbalance() #calls the function we just defined above



def getC(windms, sigmaZ,z) 
# get the concentration for a point given windspeed in meters, windms, 
#the vertical dispersion coeffecient sigmaZ, and the height in meters,z.
#using the formula described in eng2202 problem1 task handout(Ep1)( pg4/7 
#PLEASE VERIFY AND COMMENT BACK

#first broken up into 3 parts so as to end up with  c = c1 * c2 * Math.exp(cform3))
  cform1=0.79788456080286535587989211986876#(2/Math::PI)**0.5
  
#removed for efficiency
	  #cform2 = E / (sigmaZ* windms) 
	  #  $cform3div = (2*((sigmaZ)*(sigmaZ)))

  cform3 = (z*z*-1)/ $cform3div   
  c = cform1 * ($cform2 * Math.exp(cform3))
  return c  #return the computed concentration 
end

#slow implementation of the above 
#PLEASE VERIFY
def getSampleslow(x, z, metrics,sigmaZ=20)
         $cform2 = E / (sigmaZ* metrics[8])
         $cform3div = (2*((sigmaZ)*(sigmaZ)))
return getC(metrics[8], sigmaZ,z)
end

#when you run this program you can either do it from windows or command line, the command window is known as console  
def dumpInfo(met) #just that and to the console to boot. 
#make more of these to suit the format of information you want
#	 c:\dispersiondir\ruby dispersion.rb >> mydata.txt
#will capture the command output to mydata.txt
#dont forget you can get to all $variables from here as well.
print "#{met[0]} #{met[1]} windms: #{met[3]} thermal: #{met[5]} stb: #{met[7]} c@200,4: #{getSampleslow(200,4,met,getSigmaZ(4)[met[7]])}" 
end


def getSigmaZ(x)
#this gets a list of sigmaZvalues using vertical distance from source
#should probably be optimized
sigmaZ = Hash.new 
  ('A'..'F').each do |stb|
	  sigmaZ[stb]=($abcdf[stb][1] * (x ** $abcdf[stb][2]))+  $abcdf[stb][3]
  end
 # sigmaZ.each {|key, value| print key, " is ", value, "\n" }
  return sigmaZ
end

def getSample(x, z, metrics,sigmaZ=20)
#this used can be used to do some processing

#stb = metrics[7]  #windms =met[8]#sigmaZ= sigmaZ[metrics[7]]
  return getC(metrics[8], sigmaZ,z) #returns a concentration from getC bases on wind speed sigmaZ and vertical displacement
end



#GLUT is a (meta?)library for OpenGL an 3d drawing library
#basically GLUT makes drawing stuff easier and thats good.
#this draws out plots in GLUT, accepting ranges of x,z and metereological format as documented in the code  =)  
#and steps to use when traversing the timeslices
#PLEASE VERIFY
def makeframe (x,z,met,step=1)
#draw a line at x=200m
	z.each do |iz| #for each value in z as iz    
		GL.Begin(GL::POINTS); #let glut know we want to draw
		GL.Color(1, 0,1 ,1); #define the initial color for drawing
		GL.Vertex($poi[0], iz, 0,1);  #draws a point using the value defined in $poi[0] (typically a point of interest) and iz as vertical displacement in this case 
		GL.End();
	end

#get concentration for $poi 
csample = getSampleslow($poi[0],$poi[2],met,getSigmaZ($poi[2])[met[7]])
#break if csample > 0.008  #dont draw if concentration at poin was over 

#these limits are hoped to be used in depicting 
#higher concentrations with different colors rather than just
#intensity 
#ALL INPUT TO HOW TO DO SO NICELY WANTED
lowlimit = $climit#csample*csample*1
#lowlimit =  0.00003 if csample < 0.002
#(amult=100; lowlimit =  0.0004) if csample > 0.004

r,g,b,t=1,0,0,1 #default values (not used)
b = 0.002 if csample >0.002

GL.Begin(GL::POINTS);
	0.step(x.last,$step) do  |i| #for every $stepth i from 0 to the last value in x
		sigmaZ = getSigmaZ(i) #get the sigmaZ values for this x value
		sigmaZ = sigmaZ[met[7]] #we only want the sigmaZ for prevailing conditions
	
		#uses the value for sigmaZ to produce $cform2 and cform3div 
		$cform2 = E / (sigmaZ* met[8]) 
      $cform3div = (2*((sigmaZ)*(sigmaZ)))
	 
	#get a sample from this x and z=4
	#c=getSample(i,4,met,sigmaZ)
	#next  if c < lowlimit  	not used 

		0.step(z.last,$step)do |j| #for every $stepth j from 0 to the last value in z
			#get a sample at i,j
			c=getSample(i,j,met,sigmaZ)
		
			#lets c control Blue as well as concentration 
			#if c>$climit
			if c>$climit then o=1 else o=0 end
	
			GL.Color(1, 0,o*c*1000 ,c*1000);#defines the color to paint with
			GL.Vertex(i, j, 0,1);  #tells GLUT to put a dot at x,y,z,something
		end
	#drawConcPoint(i,j,c) if c>0.0001#;puts "at x,z,stb:#{i},#{j},#{met[7]} c:#{getSnapShot(i, j, met,sigmaZ)}"
	end
    
GL.End();
end


def selection(metf, ranges) #odds and ends, with metf representing the metereological datafile 
#this massages the data into suitable formats for the program 
 ranges["count"].each do #well, should do.. grab the timeslice intevals you want, for now just a oneshot loop
    met = metf.readline.split #grab a line from the metf file, and split it into a list based on space 
#the 8th (yeah the list index starts from 0, and is a clown.) 
#value is presumed to be the contrived data i gave for stability classes 
#denoted as a value from 0 to whatever
#BUT if its over 5 its reset to 5 obviously a kludge!
#need better stability class models based on windms and available data on incoming solar radiation
    met[7]="5" if met[7].to_i > 5 
    met[7].tr!('543210','ABCDEF') #translate 5 to A,  4 to B, 3 to .....etc 
#windspeed is computed and data structure is laid bare
   windkmh = met[3]
   met[8] = windkmh.to_f/3.6


    dumpInfo(met) #dump info to console/command
makeframe(ranges["x"],ranges["z"],met,10) #draw an image
  end
  
#  return met
end

def doModel() #do a model 
  ranges = Hash.new
  ranges["x"]= 0..500#distance downwind in m
  ranges["z"]= 1..300  #vertical offset
  ranges["count"]= 1..1 #should be timeslices but not used yet

selection($metf,ranges) #ranges to observe
end



#if you are reading this and you have no idea why, have a peek 
#at doModel and understand that it is called everytime T 
#is pressed or the window needs redrawing
#then keep reading upwards following functions that are called

# everything below here is glut junk





#/*  Initialize alpha blending function.
# */
def init
   GL.Enable(GL::BLEND);
   GL.BlendFunc(GL::SRC_ALPHA, GL::ONE_MINUS_SRC_ALPHA);
   GL.ShadeModel(GL::FLAT);
   GL.ClearColor(0.0, 0.0, 0.0, 0.0);
end

reshape = Proc.new {|w, h|
   GL.Viewport(0, 0,  w, h);
   GL.MatrixMode(GL::PROJECTION);
   GL.LoadIdentity();
   #if (w <= h) 
      GLU.Ortho2D(0, 500, 0.0, 300);
   #else 
     # GLU.Ortho2D(0.0, 1.0*w/h, 0.0, 1.0);
  # end
}

# /* ARGSUSED1 */
keyboard = Proc.new {|key, x, y|
   case (key)
      when 't'[0],'T'[0]
         GLUT.PostRedisplay();	
      when 27  #  /*  Escape key  */
         exit(0);
   end
}



display = Proc.new {
   GL.Clear(GL::COLOR_BUFFER_BIT);
  starttime =Time.new
   
   doModel();
   print "\t",Time.new - starttime, " model time\n"
   GL.Flush();
}

  GLUT.Init
   GLUT.InitDisplayMode(GLUT::SINGLE | GLUT::RGB);
   GLUT.InitWindowSize(500, 300);
   GLUT.CreateWindow($0);
   init();
   GLUT.ReshapeFunc(reshape);
   GLUT.KeyboardFunc(keyboard);
   GLUT.DisplayFunc(display);
   GLUT.MainLoop();