plot_snap.m 21.3 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
1
2
3
%
% Main function to produce field plots and profiles for one time step
%
4
function plot_snap(max_depl=30,max_T=0,fullV=-1,part=0,reduc=1,fin=' ',print_plot=' ') 
Lena Noack's avatar
Lena Noack committed
5
6
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D

7
8
  %addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
  addpath("C:/Arbeit/Git_Source/CHIC/Octave")
9
  %addpath("C:/Arbeit/NextCloud_FUB/CHIC/Code_Backup_March2020/Octave")
Lena Noack's avatar
Lena Noack committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

  [data,input,program_dir,fname,fpath] = read_snap(fin); % read data from snapshot
%  cd(program_dir) % return to octave folder
  if (part==1)
%    fpath
%    fname
    pname = strrep(fname,'data_snapshot','data_particles');
%    [A,input,program_dir,fpath] = read_particles;
    A = read_particles(fpath,pname);
    A(:,1:2) = A(:,1:2) * data.D / 1000.0;
  else
    A=0;
  endif  
  cd(program_dir) % return to octave folder
  [data.p_lat,data.p_rad,data.p_xx,data.p_zz] = read_mesh(fpath,input,data.nl,data.ny,data.nr); % read mesh data (position of cells)
  cd(program_dir) % return to octave folder
  data = get_data(data,input); % get fields for plotting and dimensionalize if needed
  cd(fpath); 
  mkdir "plots";
29
30
  graphics_toolkit gnuplot
  %graphics_toolkit("fltk")
31
  
32
  %size(data.depl)
33
  
Lena Noack's avatar
Lena Noack committed
34
  res=1200; %600;
Lena Noack's avatar
Lena Noack committed
35
36
37
38
39
  ar=(max(max(data.p_xx_dim))-min(min(data.p_xx_dim)))/(max(max(data.p_zz_dim))-min(min(data.p_zz_dim)));

  time = data.time/data.t_yr/1000000.0
  iter = data.iter

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
if print_plot~=' '
   name = print_plot;%z.B. 'TDHW';
   filename=strcat(name,num2str(data.iter,"%07.0f"));

   % instead of plot, print file that can be read-in and plotted in python:
   fid=fopen(sprintf('PlotData_%s.txt',filename),'w');
   fprintf(fid,'%5.10f\n',time)
   fprintf(fid,'%d %d\n',size(data.p_xx_dim))
   fprintf(fid,'%3.8f ',data.p_xx_dim)
   fprintf(fid,'\n%d %d\n',size(data.p_zz_dim))
   fprintf(fid,'%3.8f ',data.p_zz_dim)
   
   for i=1:length(name)
     if (strfind(name, "T")==i), fprintf('T'), field(:,:)=data.T(:,1,:);,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "D")==i), fprintf('D'), field(:,:)=data.depl(:,1,:);,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "H")==i), fprintf('H'), field(:,:)=data.heat(:,1,:);,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "W")==i), fprintf('W'), field(:,:)=data.w(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "v")==i), fprintf('v'), field(:,:)=data.vel(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "V")==i), fprintf('V'), field(:,:)=log10(data.visc(:,1,:)); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "a")==i), fprintf('a'), field(:,:)=data.ref_a(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "R")==i), fprintf('R'), field(:,:)=data.MP(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "c")==i), fprintf('c'), field(:,:)=data.Cp(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "k")==i), fprintf('k'), field(:,:)=data.cond(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "e")==i), fprintf('e'), field(:,:)=data.entropy(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "d")==i), fprintf('d'), field(:,:)=data.ref_r(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "M")==i), fprintf('M'), field(:,:)=data.vr(:,1,:).*data.ref_r(:,1,:)/100; ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
     if (strfind(name, "t")==i), fprintf('t'), 
       t_profile = plot_x_profile(data.T,data,input,"PrT","Temperature","southwest"," [K]");
       T_var = data.T - reshape(kron(ones(size(data.T,1),1,1),t_profile(:,2)'),size(data.T));
       field(:,:)=T_var(:,1,:); ,fprintf(fid,'\n%d %d\n',size(field)),fprintf(fid,'%3.8f ',field),endif
   endfor
   fprintf(' finished\n')

   fclose(fid)
  
else
    
  if (reduc==0) 
    % plot all relevant radial profiles (min/mean/max profiles)
79
    t_profile = plot_x_profile(data.T,data,input,"PrT","Temperature","southwest"," [K]");
80
81
82
83
    if (input.use_melt==1), t_profile = plot_x_profile(data.T,data,input,"PrTM","Temperature and Melting temp.","southwest"," [K]",log=0,melt=1);, endif
    v_profile = plot_x_profile(data.vel,data,input,"PrVel","Velocity","northeast"," [cm/yr]");
    vi_profile = plot_x_profile(data.visc,data,input,"PrVisc","Viscosity","northeast"," [log_{10} Pa\ s]",log=1);
    if(size(data.depl,1)>1), d_profile = plot_x_profile(data.depl,data,input,"PrD","Depletion","southeast","");, endif
Lena Noack's avatar
Lena Noack committed
84
85
86
87
  endif

  % plot fields
  if(max_depl==0),max_depl=max(max(data.depl));, endif
88
  min_T=0;%min(min(data.T));
Lena Noack's avatar
Lena Noack committed
89
  if(max_T==0),max_T=max(max(data.T));, endif
90
91


Lena Noack's avatar
Lena Noack committed
92
  
93
  if (reduc<=2) 
94
95
96
97
    % useful colormaps are for example hot, copper, jet (also viridis, but load extra?)
    produce_plot('T2D_',"Temperature in K",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV)
    produce_plot('Vi2D_',"Viscosity in log_{10} Pa\ s",log10(data.visc),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
    produce_plot('V2D_',"Velocity in cm/yr",data.vel,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
98
  endif
99
  
100
  if (reduc==0) 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    % data.cr_age only surface values, need to be set everywhere where there is crust
    % data.oc_vol_tot is total volume per lateral/perpendicular cell, need to get depth of crust
    nl = size(data.T,1);
    nr = size(data.T,3);
    cr_depth = zeros(nl,1);
    cr_vol = data.oc_vol_tot;
    dl = data.p_lat(2)-data.p_lat(1)
    dr = data.p_rad(2)-data.p_rad(1)
    for i=1:nl % in 3D would need do add for-loop for size(data.T,2)
      for k=nr:-1:1
        vol = 2*pi*dl*((data.p_rad(k)+dr/2)^2-(data.p_rad(k)-dr/2)^2); % vol of cell in 2D
  %      vol = 4/3*pi*dl*((data.p_rad(k)+dr/2)^3-(data.p_rad(k)-dr/2)^3); % vol of cell
        if (cr_vol(i)>0);%vol)
          cr_vol(i) = cr_vol(i)-vol;
          cr_depth(i) = cr_depth(i) + 1; % number of cells from top downwards filled as crust
        else
          break
        end
      end  
    end
    fprintf('Crust is between %d and %d cells thick\n',min(cr_depth),max(cr_depth))
    crust = zeros(size(data.T));
    for i=1:size(data.T,1) % in 3D would need do add for-loop for size(data.T,2)
      crust(i,1,nr-cr_depth(i):nr)=max(data.cr_age)-data.cr_age(i);
    end
    
    produce_plot('Cr2D_',sprintf(strcat("Crust age at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),crust,input,data,res,ar,50,'bgry',min_z=-1,max_z=-1,full=fullV)
128
129
  endif  

130
  if (reduc<=2)
131
132
133
134
135
    if(part==1)
      produce_plot('OP2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,21))
  %%    produce_plot('PD2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
  %%    produce_plot('TD2D_',"Tracer",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
    elseif(size(data.output_tracer,1)>1)
136
137
      produce_plot('O2D_',"Buoyancy (-1e10)",-data.output_tracer*1e-10,input,data,res,ar,50,'jet',min_z=0,max_z=7,full=fullV)
#      produce_plot('O2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
138
    endif
Lena Noack's avatar
Lena Noack committed
139
  endif
140
141
  
  if (reduc==0) 
142
143
144
145
146
147
148
149
    %min_z=-max(abs(min(min(data.vl))),max(max(data.vl)));
    %max_z=max(abs(min(min(data.vl))),max(max(data.vl)));
    produce_plot('Vlat2D_',"Lateral velocity in cm/yr",data.vl,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
    %min_z=-max(abs(min(min(data.vr))),max(max(data.vr)));
    %max_z=max(abs(min(min(data.vr))),max(max(data.vr)));
    produce_plot('Vrad2D_',"Radial velocity in cm/yr",data.vr,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
    if(min(min(data.strR))>0), produce_plot('StR2D_',"Strain rate in log_{10} 1/s",log10(data.strR),input,data,res,ar,50,'autumn',min_z=-1,max_z=-1,full=fullV), endif
    
Lena Noack's avatar
Lena Noack committed
150
  endif
151
  if(size(data.c,1)>1), produce_plot('C2D_',"Composition",data.c,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
152
  
Lena Noack's avatar
Lena Noack committed
153
154
155
156
157

  if(size(data.ref_r,1)>1)
    if (reduc==0) 
    produce_plot('De2D_',"Density",data.ref_r,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
    produce_plot('Mrad2D_',"Radial mass transport in kg/(m^2 yr)",data.vr.*data.ref_r/100,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
158
159
    %elseif (reduc==1)
    %produce_plot('Mradp2D_',"Total radial mass transport in kg/(m^2 yr)",abs(data.vr.*data.ref_r/100),input,data,res,ar,10,'hot',min_z=-1,max_z=-1,full=fullV)
Lena Noack's avatar
Lena Noack committed
160
161
162
    endif
  endif

163
164
  if (reduc<=2) 
    if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
165
166
  endif
  if (reduc==0) 
167
    if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
168
    if(size(data.depl,1)>1), produce_plot('D2Dc_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),-data.depl,input,data,res,ar,50,'jet',min_z=-max_depl,max_z=0,full=fullV), endif
169
170
  endif
  
Lena Noack's avatar
Lena Noack committed
171
  if (reduc==0) 
172
    if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-150,max_z=0,full=fullV), endif %input.H2O), endif
173
174
175
176
177
178
179
180
181
182
183
    if(size(data.DeltaF,1)>1), produce_plot('DF2D_',"Melt fraction",data.DeltaF,input,data,res,ar,20,'jet',min_z=-1,max_z=-1,full=fullV), endif

    if(size(data.ref_a,1)>1), produce_plot('Al2D_',"Heat expansion",data.ref_a,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.MP,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif %input.nr_ph
  %  if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=3,max_z=7,full=fullV), endif %input.nr_ph
    if(size(data.Cp,1)>1), produce_plot('Cp2D_',"Heat capacity",data.Cp,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    if(size(data.cond,1)>1), produce_plot('K2D_',"Conductivity",data.cond,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    %if(size(data.visc_diff,1)>1), produce_plot('ViDiff2D_',"Viscosity Diff in log",log10(data.visc_diff),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    %if(size(data.CO2,1)>1), produce_plot('CD2D_',"CO2 content",data.CO2,input,data,res,ar,50,'jet',min_z=0,max_z=500,full=fullV), endif %input.CO2), endif
    if(size(data.entropy,1)>1), produce_plot('En2D_',"Entropy",data.entropy,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  endif
Lena Noack's avatar
Lena Noack committed
184

185
  if (reduc==0) 
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
    % lateral temperature variations
%    data.T_var = data.T; %initialization size
%    %[tprof,tprof] = meshgrid(t_profile(:,2)); % ToDo: speed-up if t_profile matrix set up, and then just subract matrix from matrix (no loop needed)? 2D/3D - how to do that?
%
%    k=1;
%     do
%      data.T_var(:,:,k) = data.T_var(:,:,k) - t_profile(k,2); 
%      k=k+1;
%     until (k>data.nr)
%
%  %  min_z=-max(abs(min(min(data.T_var))),max(max(data.T_var)));
%  %  max_z=max(abs(min(min(data.T_var))),max(max(data.T_var)));
%    produce_plot('Tvar2D_',"Lateral temperature variation in K",data.T_var,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)%,min_z,max_z)
%
    T_var = data.T - reshape(kron(ones(size(data.T,1),1,1),t_profile(:,2)'),size(data.T));
#    size(T_var)
#    size(data.T)
#    size(kron(ones(size(data.T,1),1,1),t_profile(:,2)'))
    produce_plot('Tvar2D_',"Lateral temperature variation in K",T_var,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)%,min_z,max_z)
  endif 
Lena Noack's avatar
Lena Noack committed
206
207

  
208
  if (reduc==0) 
209
210
211
212
213
214
215
216
217
218
    if(size(data.heat,1)>1), produce_plot('H2D_',"Heat souce redistribution",data.heat,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  endif

  if (reduc==0) 
    if(size(data.U238,1)>1), produce_plot('HS_U2382D_',"U238 in ppm",data.U238,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    if(size(data.U235,1)>1), produce_plot('HS_U2352D_',"U235 in ppm",data.U235,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    if(size(data.Th232,1)>1), produce_plot('HS_Th2322D_',"Th232 in ppm",data.Th232,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
    if(size(data.K40,1)>1), produce_plot('HS_K402D_',"K40 in ppm",data.K40,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
  endif

Lena Noack's avatar
Lena Noack committed
219
220
221
  %cd(program_dir) % return to octave folder
  %colormap('default') % re-set color map
  
222
223
224
225
226
227
228
  % in case something is still open:
  fclose("all");
  %close all
  %clear all
endif %print_plot

end %function
Lena Noack's avatar
Lena Noack committed
229
230
231
232

% plot min/mean/max profiles for field x
function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dimen,log=0,melt=0)

233
234
235
236
% ATTENTION: melt profiles and pressure profile not correct, since mantle-averaged density is used
% INSTEAD: should read in profs.res and take pressure profile from there 
% (not for melt temp though, since calculated differently in convection simulationes) !!!

Lena Noack's avatar
Lena Noack committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
  nl=data.nl; % for profiles independent of periodic or reflective boundary conditions
  ny=data.ny;
  nr=data.nr;

  x_profile = zeros(nr,3);
  m_profile = zeros(nr,2);
  r = data.p_rad_dim; % radius vector

  k=1;
  do
    x_profile(k,1) = min(field(1:nl,1,k));
    x_profile(k,2) = mean(field(1:nl,1,k)); % ToDo: include dVi to get right mean value for non-uniform lateral cell distribution
    x_profile(k,3) = max(field(1:nl,1,k));
    if (melt==1)
     p = (max(r)-r(k))*input.rho*input.g/10.0^6; %km -> GPa
252
253
254
255
256
257
258
     
     if (p<12)
       DeltaTs = (0.1-input.Fe_Mantle/100.0)*(102.3 + 64.1*p - 3.62*p^2);
     else
       DeltaTs = (0.1-input.Fe_Mantle/100.0)*360.0;
     endif
     
Lena Noack's avatar
Lena Noack committed
259
     if (p<input.mz)
260
      m_profile(k,1) = input.m1s0 + input.m1s1*p + input.m1s2*p^2 + input.m1s3*p^3 + DeltaTs;
Lena Noack's avatar
Lena Noack committed
261
262
      m_profile(k,2) = input.m1l0 + input.m1l1*p + input.m1l2*p^2 + input.m1l3*p^3;
     else
263
      m_profile(k,1) = input.m2s0 + input.m2s1*p + input.m2s2*p^2 + input.m2s3*p^3 + input.m2s4*p^4 + DeltaTs;
Lena Noack's avatar
Lena Noack committed
264
265
266
267
268
269
270
      m_profile(k,2) = input.m2l0 + input.m2l1*p + input.m2l2*p^2 + input.m2l3*p^3 + input.m2l4*p^4;
     endif
    endif
    k=k+1;
  until k>nr  

  if (log==1), x_profile = log10(x_profile);, endif
271
272

  % problem: rho is mantle-averaged, hence if strongly compressible mantle, then this shifts buoyant melt radius slightly
Lena Noack's avatar
Lena Noack committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
  
  f = figure('visible','off');
  
  h=plot(x_profile(:,1),r,"b");
  hold on
  h2=plot(x_profile(:,2),r,"k");
  h3=plot(x_profile(:,3),r,"r");
  if (melt==1)
    h4=plot(m_profile(:,1),r,":g");
    h5=plot(m_profile(:,2),r,"--g");
  endif
  set(h,'LineWidth',5)
  set(h2,'LineWidth',10)
  set(h3,'LineWidth',5)
  if (melt==1)
    set(h4,'LineWidth',5)
    set(h5,'LineWidth',5)
  endif
291
292
293
294
295
296
297
298
  
  if (melt==1)
    radius_melt_buoyant = max(r)-input.md*10.0^6/(input.rho*input.g);
    plot([min(x_profile(:,1)),max(x_profile(:,3))],[radius_melt_buoyant,radius_melt_buoyant],'k:')
    
    axis([min(x_profile(:,1)),3000,max(r)-3*input.md*10.0^6/(input.rho*input.g),max(r)])
  endif
  hold off
Lena Noack's avatar
Lena Noack committed
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
%  legend(sprintf(strcat(name,' min')),sprintf(strcat(name," mean")),sprintf(strcat(name," max")), "location", orient);
%  legend("boxoff");
  
  FN = findall(f,'-property','FontName');
  set(FN,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf');
  FS = findall(f,'-property','FontSize');
  set(FS,'FontSize',20);

%  dimen="";
%  if (data.D>1.0), dimen=" [K]";, endif
  dimenR="";
  if (data.D>1.0), dimenR=" [km]";, endif

  xlabel (sprintf(strcat(title_plot,dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',24);
  ylabel (sprintf(strcat("Radius",dimenR)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',24);
  title (sprintf(strcat(title_plot," profile")),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Bold.ttf','FontSize',30);
  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);                                                                          

%  graph=strcat("-S", num2str((ar+0.4)*res),",", num2str(res));
318
%  print(filename,'-deps','-S800,720');%graph); 
Lena Noack's avatar
Lena Noack committed
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
  print(filename,'-dpng','-S800,720');%graph); 

  clear f;
  
end

% help function for plot_snap to produce png plots for selected fields
function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,min_z=-1,max_z=-1,full=-1,partX=0,partY=0,partV=0)

%  min_z
%  max_z

  if (min_z==-1)
    min_z = min(min(field));
  endif
  if (max_z==-1)
    max_z = max(max(field));
  endif

338
  font_title=round(20*res/600);%(30*res/600);
Lena Noack's avatar
Lena Noack committed
339
340
341
342
  font_axes_title=round(24*res/600);
  font_axes=round(20*res/600);
  
  f = figure('visible','off');
343

Lena Noack's avatar
Lena Noack committed
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
  if (data.ny>1)
    ny_use=data.ny/2;
  else
    ny_use=1;
  endif
  dummy(:,:)=field(:,ny_use,:);
%  contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour); % works both for spherical and cartesian plots

  %[size(data.p_xx_dim),  size(data.p_zz_dim),  size(dummy), size(field)]
%  size(partX,1)
  if (size(partX,1)>1) 
%    fprintf("Particle plot")
%    size(partX,1)
    N=length(partX);
    colorMap = zeros(N, 3);
%    for k = 1:N
%     if (input.set_tracers==1)
%     elseif (input.set_tracers==2)
%%       colorMap(k,:) = (partV(k)-min(partV))/(max(partV)-min(partV))*[1,0,0]+(1-(partV(k)-min(partV))/(max(partV)-min(partV)))*[0,0,1]; % shifts data points between red and blue depending on value
%       colorMap(k,:) = partV(k)/max(partV)*[1,0,0]+(1-partV(k)/max(partV))*[0,0,1]; % shifts data points between red and blue depending on value
%     endif
%    end
    
    if (input.set_tracers==1)
     %colorMap(partV<0.25)
     for k = 1:N
      if (partV(k)<0.25)
        colorMap(k,:) = [1,0,0];
      elseif (partV(k)>0.75)
        colorMap(k,:) = [1,1,1];
      else
        colorMap(k,:) = [0,0,1];
      endif
     end        
    elseif (input.set_tracers==2)
      if (max(partV) != min(partV))
        colorMap(:,1) = (partV-min(partV))/(max(partV)-min(partV));
        colorMap(:,2) = 1 - (partV-min(partV))/(max(partV)-min(partV));
      else
        colorMap(:,1) = partV/max(partV);
        colorMap(:,3) = 1 - partV/max(partV);
      endif
    end
    
%    dummy(:,:)=0.4;
%    contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
%    hold on
    %scatter(partX,partY,100,colorMap,'filled')
    scatter(partX,partY,100,partV,'filled')
    %partX
    %plot(partX,partY)
%    hold off
  else  
%    contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
    [C,h] = contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
399
400

%    imagesc(data.p_lat,data.p_rad,dummy);
Lena Noack's avatar
Lena Noack committed
401
    set(h,'LineColor','none')
402
    set(gca,'visible','off')
Lena Noack's avatar
Lena Noack committed
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
  endif

%  surface(data.p_xx_dim,data.p_zz_dim,dummy,'EdgeColor','none'); % works both for spherical and cartesian plots

%  surfc(data.p_xx_dim,data.p_zz_dim,dummy,'EdgeColor','none'); % works both for spherical and cartesian plots
%  shading interp;
%  view(0,90)


% levelplot
  if (strcmp(cmap,'bgry'))
    mycolormap % call mycolormap function to get colormap (mymap) black-blue-green-red-yellow
    colormap(mymap)
  else
    colormap(cmap)
  endif
419
420
421
422
423
%    min_z
%    max_z
  if (min_z<max_z)
    caxis([min_z max_z]);
  endif
Lena Noack's avatar
Lena Noack committed
424
425
426
427
428
429
430
431
432
433
434
435
  colorbar

  FN = findall(f,'-property','FontName');
  set(FN,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf');%dejavu/DejaVuSans.ttf');
  FS = findall(f,'-property','FontSize');
  set(FS,'FontSize',font_axes);

%  axis ("xy", "square");  
%  axis ([0, ar, 1, 2]);
  axis ([min(min(data.p_xx_dim)), max(max(data.p_xx_dim)), min(min(data.p_zz_dim)), max(max(data.p_zz_dim))]); % min_z max_z dazu?
  dimen="";
  if (data.D>1.0), dimen=" [km]";, endif
436
437
%  xlabel (sprintf(strcat("Lateral direction",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
%  ylabel (sprintf(strcat("Radius",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
Lena Noack's avatar
Lena Noack committed
438
439
440
441
442
443
  title (title_plot,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Bold.ttf','FontSize',font_title);%dejavu/DejaVuSans-Bold.ttf','FontSize',font_title);

  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);
%  filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.eps'),f);

  if (full==-1)
444
  %  fprintf("Regional Sphere")
Lena Noack's avatar
Lena Noack committed
445
446
    graph=strcat("-S", num2str((ar+0.6)*res),",", num2str(res)); % suitable for regional sphere
  else
447
  %  fprintf("Full Sphere")
Lena Noack's avatar
Lena Noack committed
448
449
    graph=strcat("-S", num2str((ar+0.37)*res),",", num2str(res)); % suitable for full sphere
  endif
450
  %print(filename,'-deps',graph);   
Lena Noack's avatar
Lena Noack committed
451
452
453
454

  print(filename,'-dpng',graph);

  clear f;
Lena Noack's avatar
Lena Noack committed
455
end