lib/gs2crmod/calculations.rb in gs2crmod-0.12.15 vs lib/gs2crmod/calculations.rb in gs2crmod-0.12.16

@@ -238,13 +238,10 @@
 # 	
 # 	ep t_vec.dup.delete_if{|el| not (el - 71.3).abs < 0.5}
 # 	
 # 	exit
 	
-	
-	
-	
 	return
 	
 	# Calculate a series of time averaged segments
 	pieces = hflux.pieces(20) # split into 20 pieces
 	avgs = GSL::Vector.alloc(pieces.map{|vec| vec.sum/vec.size})
@@ -270,21 +267,22 @@
 
 alias :csti :calculate_saturation_time_index
 
 def calculate_frequencies
   @real_frequencies = FloatHash.new
+  @frequency_at_ky_at_kx ||= FloatHash.new
   omega_avg_narray = netcdf_file.var("omega_average").get('start' => [0, 0, 0, -1], 
                                                           'end' => [0, -1, -1, -1])
   omega_avg_narray.reshape!(*omega_avg_narray.shape.slice(1..2))
 
   if @grid_option == 'single'
     @frequency_at_ky_at_kx = omega_avg_narray[0]
   else
     list(:ky).values.sort.each_with_index do |kyv, i|
       @frequency_at_ky_at_kx[kyv] = FloatHash.new
       list(:kx).values.sort.each_with_index do |kxv, j|	
-        @frequency_at_ky_at_kx[kyv][kxv] = omega_avg_narray[i, j]
+        @frequency_at_ky_at_kx[kyv][kxv] = omega_avg_narray[j, i]
       end
       write_results
     end
   end
 end
@@ -369,28 +367,27 @@
 end
 
 alias :cgrf :calculate_growth_rates_and_frequencies
 
 def calculate_growth_rate(vector, options={})
-	raise "This vector should be positive definite" if vector.min < 0.0
-	offset = 0
-	length = vector.length
+  raise "This vector should be positive definite" if vector.min < 0.0
+  offset = 0
+  length = vector.length
   while vector[offset] == 0.0
     offset+=1
     return 0.0 if offset == vector.length
   end
-	growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length - offset)))[1]
-	divisor = 1
-	while (growth_rate.to_s == "NaN")
-      #This corrects the growth rate if phi has grown all the way to NaN during the simulation
-      divisor *= 2
-      length = (vector.size.to_f / divisor.to_f).floor
-      return "NaN" if length <= offset + 1
-      growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length-offset)))[1]
-	end	
-    p 'HELLLOOASJDKFHASDF'
-	growth_rate
+  growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length - offset)))[1]
+  divisor = 1
+  while (growth_rate.to_s == "NaN")
+    #This corrects the growth rate if phi has grown all the way to NaN during the simulation
+    divisor *= 2
+    length = (vector.size.to_f / divisor.to_f).floor
+    return "NaN" if length <= offset + 1
+    growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length-offset)))[1]
+  end	
+  growth_rate
 end
 
 # Not needed for GS2 built after 16/06/2010
 def corrected_mom_flux_stav
 	par_mom_flux_stav - perp_mom_flux_stav
@@ -681,20 +678,20 @@
 	end
 end
 
 def max_trans_phi
 	phivec = gsl_vector('phi2tot_over_time')
-	offset = 30
+	#offset = 30
 	phivec.subvector(20, phivec.size - 20).max
 end
 
 def max_es_heat_amp(species_index)
 	@transient_es_heat_flux_amplification_at_species_at_ky[species_index-1].values.max
 end
 
 def calculate_spectral_checks
-  kx = gsl_vector('kx')
+  #kx = gsl_vector('kx')
   ky = gsl_vector('ky')
   ky_spec = gsl_vector('spectrum_over_ky')
   kx_spec = gsl_vector('spectrum_over_kx')
   kpar_spec = gsl_vector('spectrum_over_kpar', ky_index: ky_spec.max_index + 1, kx_index: 1)
   
@@ -750,10 +747,17 @@
 	return @spectrum_check.min >= min
 end
 alias :csc :calculate_spectral_checks
 
 def calculate_prandtl_number
+  
+  @prandtl_number = {}
+  @nspec.times do |i|
+    @prandtl_number[i+1] = nil
+  end
+  write_results
+
   if @g_exb == 0
     eputs 'g_exb = 0 therefore Prandtl number is undefined.'
     return nil
   elsif @nonlinear_mode=="off"
     eputs 'Prandtl number only makes sense for a nonlinear run.'
@@ -761,12 +765,10 @@
   elsif @local_eq.fortran_false?
     eputs 'Prandtl number currently only calculated for Miller equilibrium.'
     return nil
   end
 
-  @prandtl_number = {}
-
   @nspec.times do |i|
     species_index = i + 1
     @prandtl_number[species_index] = - (@rhoc/@qinp/@rmaj**2) * 
                                        (eval("@tprim_#{species_index}")/@g_exb) * 
                                        (@es_mom_flux_stav[species_index]/
@@ -775,6 +777,5 @@
   write_results
 end
 
 end
 end
-