Changeset 38
- Timestamp:
- 01/16/08 01:52:08 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/procs/plt_map.pro
r34 r38 938 938 939 939 'yfx': BEGIN 940 ; 941 ; Scatter plot y=f(x) 942 ; ------------------- 943 ; 944 t = tmask 945 boxx = '' 946 boxy = '' 947 add_txt = '' 948 ; arrange/average data if needed 949 ; space and time plots fld 950 IF strpos(cmd.plt, 'xy_') NE -1 THEN BEGIN 951 mask_z, fld, cmd, boite_plt, dimplot, boxx 952 boxx = ' ['+boxx+']' 953 print, ' Averaging (space) '+cmd.var+' in'+boxx 954 fld = checkfield(fld, 'plt', type = 'xy', boite = boite_plt) 955 ENDIF 956 IF datyp.hotyp NE '-' THEN BEGIN 957 mask_z, fld, cmd, boite_plt1d, dimplot, boxx 958 boxx = ' ['+boxx+']' 959 vargrid = vargrid1 960 IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d[0:3], vargrid 961 domdef, boite_plt1d[0:3] 962 print, ' Averaging (time serie plot) cmd.var '+cmd.var+' in'+boxx 963 IF debug_w THEN print, 'fld size', size(fld) 964 fld = checkfield(fld, 'plt1d', type = datyp.hotyp, boite = boite_plt1d) 965 IF cmd.trend GT 0 THEN BEGIN 966 fld_flag = 1 967 fld = trends(fld, cmd.trend, datyp.hotyp) 968 add_txt = trd_txt 969 ENDIF 970 ENDIF 971 ; space and time plots fld2 972 IF sw_diffg EQ 1 THEN BEGIN 973 jptb = jpt 974 IF debug_w THEN print, ' calling def_grid, cmd2 dans plt_map ' 975 def_grid, cmd2 976 jpt = jptb 977 ENDIF 978 IF strpos(cmd2.plt, 'xy_') NE -1 THEN BEGIN 979 mask_z, fld2, cmd2, boite_plt2, dimplot, boxy 980 boxy = ' ['+boxy+']' 981 print, ' Averaging (space) '+cmd2.var+' in'+boxy 982 fld2 = checkfield(fld2, 'plt', type = 'xy', boite = boite_plt2) 983 ENDIF 984 IF datyp2.hotyp NE '-' THEN BEGIN 985 run_stddev = 0 986 IF strpos(cmd2.plt, '@r') GE 1 THEN BEGIN 987 run_stddev = strmid(cmd2.plt, strpos(cmd2.plt, '@r')+2, strlen(cmd2.plt)) 988 cmd2.plt = strmid(cmd2.plt, 0, strpos(cmd2.plt, '@r')) 989 ENDIF 990 mask_z, fld2, cmd2, boite_plt1d2, dimplot, boxy 991 boxy = ' ['+boxy+']' 992 vargrid = vargrid2 993 IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d2[0:3], vargrid 994 domdef, boite_plt1d2[0:3] 995 print, ' Averaging (time serie plot) cmd2.var '+cmd2.var+' in'+boxy 996 IF debug_w THEN print, ' fld2 size', size(fld2) 997 fld2 = checkfield(fld2, 'plt1d', type = datyp2.hotyp, direc = 'xy', boite = boite_plt1d2, /timearray) 998 IF cmd2.trend GT 0 THEN BEGIN 999 fld_flag = 2 1000 fld2 = trends(fld2, cmd2.trend, datyp2.hotyp) 1001 add_txt = trd_txt 1002 ENDIF 1003 ENDIF 1004 IF sw_diffg EQ 1 THEN BEGIN 1005 sw_diffg = 0 1006 jptb = jpt 1007 IF debug_w THEN print, ' calling def_grid, cmd dans plt_map ' 1008 def_grid, cmd 1009 jpt = jptb 1010 ENDIF 1011 coeff = linfit(fld2, fld, CHISQ = linerr, PROB = proberr, SIGMA = sigmaerr) 1012 print, ' Linfit coef (a+bx) = ', coeff(0), coeff(1), ' errbar = ',sigmaerr, ' (proba = ',proberr, ')' 1013 1014 ; plot domain 1015 boxyfx = def_box(cmd.plt, dimplot, legz, time_stride) 1016 ; define line color and type 1017 overc = overlay_type(iover, dimplot) 1018 vardate = date_txt 1019 varname = varlegend 1020 varunit = ' '+cmd.var+boxx+' =f('+cmd.var2+boxy+') '+add_txt 1021 pltcmdstd = 'pltsc,fld,fld2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+overc+com_strplt 1022 1023 IF hotyp EQ 't' THEN BEGIN ; time scatter plot 1024 IF mean_sc_only EQ 0 OR mean_sc_only EQ 4 THEN BEGIN 1025 ; number of time colors 1026 IF strpos(symbol_families, 'x') NE -1 THEN BEGIN 1027 coding = str_sep(symbol_families, 'x') 1028 ncolors = long(coding(0)) 1029 ntimes = long(coding(1)) 1030 ENDIF ELSE BEGIN 1031 ncolors = long(symbol_families) 1032 ntimes = 1 1033 ENDELSE 1034 IF ncolors GE 2 THEN BEGIN ; time color coding 1035 ic = 0 1036 jpto = jpt 1037 jpt = jpt/ncolors 1038 lincoef = fltarr(ncolors) 1039 linerro = fltarr(ncolors) 1040 linprob = fltarr(ncolors) 1041 WHILE ic LE ncolors-1 DO BEGIN 1042 idx0 = (floor(findgen(jpto/(ncolors*ntimes)))*ncolors*ntimes)+ic*(ntimes) 1043 idx = idx0 1044 jl = 1 1045 WHILE jl LE ntimes-1 DO BEGIN 1046 idx = [idx, idx0+jl] 1047 jl = jl+1 1048 ENDWHILE 1049 fldp = fld(idx) 1050 fldp2 = fld2(idx) 1051 coeff = linfit(fldp2(where (fldp2(*) GT 0.)), fldp(where (fldp2(*) GT 0.)), CHISQ = errlin, PROB = prblin) 1052 linerro(ic) = errlin 1053 linprob(ic) = prblin 1054 lincoef(ic) = coeff(1) 1055 ; print, ' Period '+strtrim(string(ic), 2)+' Linear fit slope, chisq, prob =', coeff(1)*1000., errlin*1000., prblin 1056 IF mean_sc_only EQ 0 THEN BEGIN 1057 pltcmd = 'pltsc,fldp,fldp2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d='+string(ic)+',COLOR='+string(symbol_color(ic))+', STY1D='+string(symbol_style(ic)-1)+', SYMSIZE='+string(symbol_size) 1058 printf, nulhis, strcompress(pltcmd) 1059 IF debug_w THEN print, ' ',pltcmd 1060 res = execute(pltcmd) 1061 ENDIF 1062 IF mean_sc_only EQ 4 AND ic EQ ncolors-1 THEN BEGIN 1063 print, ' Linfit slope + error (June-Nov)=', mean(lincoef(5:10)), mean(linerro(5:10)) 1064 print, ' Linfit slope + error (Dec-May)=', mean([lincoef(0:4), lincoef(11)]), mean([linerro(0:4), linerro(11)]) 1065 vardate = 'toto' 1066 varname = 'Coupling strength '+cmd.date1+' '+cmd.spec 1067 jpt = ncolors 1068 time = lindgen(12)*30*1+ $ 1069 julday(1,1,01, ndayspm = calendar_type) 1070 hotyp = 't' 1071 minmax = ',0,25' 1072 pltcmd = 'pltt,lincoef*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=0,COLOR=1, thick=4, STY1D=0' 1073 1074 printf, nulhis, strcompress(pltcmd) 1075 IF debug_w THEN print, ' ', pltcmd 1076 res = execute(pltcmd) 1077 pltcmd = 'pltt,(lincoef-linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0' 1078 printf, nulhis, strcompress(pltcmd) 1079 IF debug_w THEN print, ' ',pltcmd 1080 res = execute(pltcmd) 1081 pltcmd = 'pltt,(lincoef+linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0' 1082 printf, nulhis, strcompress(pltcmd) 1083 IF debug_w THEN print, ' ',pltcmd 1084 res = execute(pltcmd) 1085 ENDIF 1086 ic = ic + 1 1087 ENDWHILE 1088 ENDIF ELSE BEGIN ; one color 1089 IF debug_w THEN print, ' ',pltcmdstd 1090 res = execute(pltcmdstd(0)) 1091 ENDELSE 1092 ENDIF 1093 ENDIF ELSE BEGIN ; standard scatter plot 1094 IF debug_w THEN print, ' ',pltcmdstd 1095 res = execute(pltcmdstd(0)) 1096 ENDELSE 1097 1098 ; add mean SC in 4x3 plots 1099 1100 IF symbol_families EQ '4x3' THEN BEGIN 1101 1102 IF cmd.trend EQ '412' THEN BEGIN 1103 fldrem_t1 = fldrem_t1-mean(fldrem_t1) 1104 fldrem_t2 = fldrem_t2-mean(fldrem_t2) 1105 ENDIF ELSE BEGIN 1106 running = 12L 1107 lenght = (size(fld))[1] 1108 fldrem_t1 = fltarr(running) 1109 fldrem_t2 = fltarr(running) 1110 FOR t1 = 0, running-1 DO BEGIN 1111 fldrem_t1(t1) = mean(fld(long(findgen(lenght/running)*running+t1))) 1112 fldrem_t2(t1) = mean(fld2(long(findgen(lenght/running)*running+t1))) 1113 ENDFOR 1114 ENDELSE 1115 1116 print, ' Seasonal cycle var/stddev fld1 ', (moment(fldrem_t1))[1], sqrt((moment(fldrem_t1))[1]) 1117 print, ' Seasonal cycle var/stddev fld2 ', (moment(fldrem_t2))[1], sqrt((moment(fldrem_t2))[1]) 1118 1119 fldrem_t1 = [fldrem_t1, fldrem_t1(0)] 1120 fldrem_t2 = [fldrem_t2, fldrem_t2(0)] 1121 sw_ov1d = mean_sc_only 1122 IF mean_sc_only EQ 2 AND cmd.trend EQ '412' THEN BEGIN 1123 ; SC of std dev 1124 stat = spectra(fld, jpt/12, 20, 15, 0) 1125 stat2 = spectra(fld2, jpt/12, 20, 15, 0) 1126 print, ' Stddev of monthly stddevs fld', sqrt((moment(stat.sc_std-(moment(stat.sc_std))[0]))[1]) 1127 print, ' Stddev of monthly stddevs fld2', sqrt((moment(stat2.sc_std-(moment(stat2.sc_std))[0]))[1]) 1128 stdsc = [stat.sc_std, stat.sc_std(0)] 1129 stdsc2 = [stat2.sc_std, stat2.sc_std(0)] 1130 pltcmd = 'pltsc,stdsc,stdsc2,0. ,maxc,0.,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=1-min(1,sw_ov1d), STY1D=-1, THICKN=3, SYMSIZE='+string(symbol_size)+', FRACTION=fraction' 1131 printf, nulhis, strcompress(pltcmd) 1132 IF debug_w THEN print, pltcmd 1133 res = execute(pltcmd(0)) 1134 pltcmd = 'xyouts,stdsc2-(maxc2)/(.3*win(0)*win(1)),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' 1135 printf, nulhis, strcompress(pltcmd) 1136 IF debug_w THEN print, pltcmd 1137 res = execute(pltcmd(0)) 1138 ENDIF ELSE BEGIN 1139 ov1d_val = 0 1140 IF mean_sc_only EQ 0 THEN ov1d_val = 1 1141 pltcmd = 'pltsc,fldrem_t1,fldrem_t2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=ov1d_val, STY1D=-1, THICKN=5, SYMSIZE='+string(symbol_size)+', FRACTION=fraction' 1142 printf, nulhis, strcompress(pltcmd) 1143 IF debug_w THEN print, pltcmd 1144 res = execute(pltcmd(0)) 1145 IF mean_sc_only EQ 1 THEN BEGIN 1146 ; add month info win(0)*win(1) 1147 pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5*win(0)*win(1)),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' 1148 printf, nulhis, strcompress(pltcmd) 1149 IF debug_w THEN print, pltcmd 1150 res = execute(pltcmd(0)) 1151 ENDIF 1152 ENDELSE 1153 ENDIF 1154 1155 ; overplot sigma contours in T-S plane 1156 IF (cmd.var EQ 'sosstsst' and cmd.var2 EQ 'sosaline') OR (cmd.var EQ 'votemper' and cmd.var2 EQ 'vosaline') THEN BEGIN 1157 IF noerase EQ 0 THEN BEGIN 1158 min_t = minc-2.0 1159 max_t = maxc+2.0 1160 min_s = minc2-2.0 1161 max_s = maxc2+2.0 1162 delta_t = 0.5 1163 delta_s = 0.05 1164 t_vec = findgen((max_t-min_t)/delta_t+1)*delta_t+min_t 1165 s_vec = findgen((max_s-min_s)/delta_s+1)*delta_s+min_s 1166 t_ = t ; since t is already defined 1167 t = t_vec 1168 FOR i = 1, (max_s-min_s)/delta_s DO t = [[t], [t_vec]] 1169 s = s_vec 1170 FOR i = 1, (max_t-min_t)/delta_t DO s = [[s], [s_vec]] 1171 s = transpose(s) 1172 sr=sqrt(abs(s)) 1173 r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $ 1174 -9.095290E-3)*t+6.793952E-2)*t+999.842594 1175 r2=(((5.3875E-9*t-8.2467E-7)*t+7.6438E-5)*t-4.0899E-3)*t+8.24493E-1 1176 r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3 1177 rhopn = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) -1000.0 1178 min_sig = round(min(rhopn)+0.5)-1. 1179 max_sig = round(max(rhopn)-0.5)+1. 1180 delta_sig = 1. ; define contour interval here! 1181 levels = findgen((max_sig-min_sig)/delta_sig+1)*delta_sig+min_sig 1182 labels = levels 1183 labels(*) = 1 1184 contour, rhopn, s, t, /overplot, levels = levels, c_labels = labels, c_linestyle = 0 1185 t = t_ 1186 ENDIF 1187 ENDIF 1188 1189 printf, nulhis, 'boite_yfx=',boxyfx 1190 printf, nulhis, strcompress(pltcmd) 1191 1192 tmask = t 1193 ; force read grid for next window 1194 cmd_prev.grid = 'toto' 1195 940 @yfx 1196 941 END 1197 942 'spec': BEGIN
Note: See TracChangeset
for help on using the changeset viewer.