#!/usr/bin/ruby # -*- coding: utf-8 -*- # =begin # == 概要 # # 横軸に惑星半径, 縦軸に温度 をとった図を作成する # # == 更新履歴 # * 2014-12-25 増田 和孝 半径計算用に改造 # * 2013-11-15 石渡 地球放射計算用に変更 # * 2011/08/12 納多 哲史 Omega = 0.799 (0.8 で南北非対称が出た場合) を追加 # * 2011/04/27 納多 哲史 新規作成 # # =end require "numru/ggraph" require 'fileutils' include NumRu LatentHeat = 2.5e6 # x 軸をログプロットする場合 #axis_x_log = true axis_x_log = false #read_meant_data = false DATA_HOME_LOCAL = '/GFD_Dennou_Work11/masuda/' # config.yml で設定されているディレクトリ TMPFILE = "dcl.ps" TIME_START = 725 TIME_END = 1095 VAL_MAX = 200 VAL_MIN = 170 WSN = 2 MARK_SIZE = 0.01 #YNAME = 'Integrated Heat Flux' #TITLE = 'Heat Budget' TITLE = '' #1: 点 #2: 十字 #3: アスタリスク #4: 白抜き丸 #5: バツ #6: 白抜き正方形 #7: 白抜き上向き三角 #8: 白抜き菱形 #9: 白抜き星 #10: 黒塗り丸 #11: 黒塗り正方形 #12: 黒塗り上向き三角 #13: 黒塗り左向き三角 #14: 黒塗り下向き三角 #15: 黒塗り右向き三角 #16: 黒塗り星 #17: 黒塗り右向き旗 mark_type = Hash.new mark_type = 10 rplanets=[0.5,1.0,1.25,1.46,2.0,4.0,6.0,8.0,12.0,16.0] na_rplanets = NArray.to_na(rplanets) nrplanet = rplanets.size X_MIN = rplanets.min - 0.5 X_MAX = rplanets.max + 0.5 # 描画設定 #p i # if i == 0 then DCL.uzfact(0.6) # end night = NArray.sfloat(nrplanet) ## 各 rplanet で回す. for i in 0..nrplanet-1 rplanet = na_rplanets[i] p "Temp" + ' ' + rplanet.to_s ### 読み込み dir = 'E_T42L26_sync_R' + rplanet.to_s file = 'Temp.nc' path = File.join(DATA_HOME_LOCAL, dir, file) gphys = GPhys::IO.open(path, 'Temp') na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_sig = GPhys::IO.open(path,'sig').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val na_sig_weight = GPhys::IO.open(path, 'sig_weight').val nlon = na_lon.size nlat = na_lat.size nsig = na_sig.size ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>TIME_START..TIME_END).mean('time') #### 空間平均 ##### 緯度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum gphys = gphys.sum('lat') gphys = gphys * na_sig_weight.reshape(1, nsig) / na_sig_weight.sum gphys = gphys.sum('sig') night[i] = gphys.cut('lon'=>180..360).mean('lon').val # gphys = (gphys*na_lat_weight).sum(1)/na_lat_weight.sum # night[i] = gphys.val p night[i] ## check # p all - (west + east)/2.0 end ## gphys オブジェクトの生成 na_rplanets_4axis = na_rplanets na_rplanets_4axis = NMath::log10( na_rplanets ) if axis_x_log # (2013-6-14 石渡) 本当は, rplanet 軸の long_name は Rplanet/Rplanet_E # にしたかったのだけど, なぜかうまくいかない. # しかたなく, 今は Rplanet/RplanetE にしてある. tmp_long_name = "Rplanet/RplanetE" tmp_long_name = 'log10(' + tmp_long_name + ')' if axis_x_log va_rplanet = VArray.new( na_rplanets_4axis, {"long_name"=>tmp_long_name}, "rplanets" ) axis_rplanet = Axis.new.set_pos(va_rplanet) # testing gp = Hash.new long_name = 'Temperature' #units = 'W m-2' units = 'K' p long_name p units p night data = VArray.new( night, {"long_name"=>long_name, "units"=>units}, "Temp" ) gp = GPhys.new( Grid.new(axis_rplanet), data ) # 解像度を変えた実験(T170) t85n = NArray.sfloat(nrplanet) # rplanet = na_rplanets[nrplanet-1] # p "SurfTemp" + ' ' + rplanet.to_s # # 8.0 R dir = 'E_T85L26_sync_R8.0' file = 'SurfTemp.nc' path = File.join(DATA_HOME_LOCAL, dir, file) gphys = GPhys::IO.open(path, 'SurfTemp') na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>1095..1195).mean('time') #### 空間平均 ##### 緯度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum t85n[7] = gphys.cut('lon'=>180..360).mean('lon').sum('lat').val # gphys = (gphys*na_lat_weight).sum(1)/na_lat_weight.sum ## night[i] = gphys.val # p night[nrplanet-1] # ## check # # p all - (west + east)/2.0 ### 読み込み # 16.0R dir = 'E_T85L26_sync_R16.0' file = 'SurfTemp.nc' path = File.join(DATA_HOME_LOCAL, dir, file) gphys = GPhys::IO.open(path, 'SurfTemp') na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>330..660).mean('time') #### 空間平均 ##### 緯度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum t85n[nrplanet-1] = gphys.cut('lon'=>180..360).mean('lon').sum('lat').val # gphys = (gphys*na_lat_weight).sum(1)/na_lat_weight.sum ## night[i] = gphys.val # p night[nrplanet-1] # ## check # # p all - (west + east)/2.0 p t85n data85 = VArray.new(t85n,{"long_name"=>long_name,"units"=>units},"SurfTemp(T85)") gp2 = GPhys.new(Grid.new(axis_rplanet),data85) # 解像度を変えた実験(T170) t170n = NArray.sfloat(nrplanet) # rplanet = na_rplanets[nrplanet-1] # p "SurfTemp" + ' ' + rplanet.to_s # # 8.0 R dir = 'E_T170L26_sync_R8.0' file = 'SurfTemp.nc' path = File.join(DATA_HOME_LOCAL, dir, file) gphys = GPhys::IO.open(path, 'SurfTemp') na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>1155..1195).mean('time') #### 空間平均 ##### 緯度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum t170n[7] = gphys.cut('lon'=>180..360).mean('lon').sum('lat').val # gphys = (gphys*na_lat_weight).sum(1)/na_lat_weight.sum ## night[i] = gphys.val # p night[nrplanet-1] # ## check # # p all - (west + east)/2.0 data170 = VArray.new(t170n,{"long_name"=>long_name,"units"=>units},"SurfTemp(T170)") gp3 = GPhys.new(Grid.new(axis_rplanet),data170) ## 描画 # 描画準備 # testing DCL.gropn(WSN) GGraph.set_fig('window'=>[X_MIN,X_MAX,VAL_MIN,VAL_MAX],'viewport'=>[0.3,0.6,0.3,0.5]) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', false) # 初期化? # DCL.grfrm #GGraph.set_axes('xtickint'=>xlabel_int) #GGraph.set_axes('xlabelint'=>xlabel_int) #GGraph.set_axes('ytickint'=>ylabel_int) #GGraph.set_axes('ylabelint'=>ylabel_int) #GGraph.line( 'max'=>val_max, 'min'=>val_min ) # オプションのデフォルト値の設定 opt = { 'max'=>VAL_MAX, 'min'=>VAL_MIN, 'size'=>MARK_SIZE, 'title'=>TITLE,'type'=>mark_type,'legend'=>true } GGraph.mark(gp, true,opt) # eps ファイルが存在する場合は終了 # if File.exist?(epsfn) then # puts("MESSAGE: #{epsfn} have already existed. Skipped.") # next # end # testing # DCL.gropn(WSN) # p draw_type + '_' + area psfn = 'RPlanet_Temp-nightmean.ps' epsfn = 'RPlanet_Temp-nightmean.eps' # ps ファイルの処理 # remove extra space of ps file # if WSN == 2 then # if File.exist?(TMPFILE) then # debug DCL.grcls # DCL.grfrm # DCL.grfig FileUtils.mv(TMPFILE, psfn) puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrmcm #{psfn} | dclpsrot | dclpsmargin > #{epsfn}` File.delete(psfn) # testing exit # end # end # debug # DCL.grcls exit(0)