dispersion

#Dispersion.rb implements the gaussian dispersion model using the ruby language#Bugs,Comments or Questions unomi@techemail.comrequire "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 variablesE = 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 number200.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 ondef 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 abovedef 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 VERIFYdef 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]])}" enddef getSigmaZ(x)#this gets a list of sigmaZvalues using vertical distance from source#should probably be optimizedsigmaZ = 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 sigmaZenddef 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 displacementend#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 VERIFYdef 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 WANTEDlowlimit = $climit#csample*csample*1#lowlimit =  0.00003 if csample < 0.002#(amult=100; lowlimit =  0.0004) if csample > 0.004r,g,b,t=1,0,0,1 #default values (not used)b = 0.002 if csample >0.002GL.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();enddef 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/commandmakeframe(ranges["x"],ranges["z"],met,10) #draw an image  end  #  return metenddef 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 yetselection($metf,ranges) #ranges to observeend#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);endreshape = 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();