Changeset 38 for trunk/procs/plt_map.pro


Ignore:
Timestamp:
01/16/08 01:52:08 (16 years ago)
Author:
ericg
Message:

Extracted yfx.pro from plt_map.pro as it was getting too long

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/procs/plt_map.pro

    r34 r38  
    938938 
    939939            '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 
    1196941            END 
    1197942            'spec': BEGIN 
Note: See TracChangeset for help on using the changeset viewer.