[cig-commits] r8975 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Wed Jan 2 00:31:29 PST 2008
Author: becker
Date: 2008-01-02 00:31:28 -0800 (Wed, 02 Jan 2008)
New Revision: 8975
Added:
mc/1D/hc/trunk/manipulate_hc.py
mc/1D/hc/trunk/pdens.py
mc/1D/hc/trunk/pvisc.py
Modified:
mc/1D/hc/trunk/README.TXT
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_init.c
Log:
Modified the way that the surface boundary conditions are handled on the command line.
Added python scripts to plot and interactively modify viscosity and density scaling profiles.
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/README.TXT 2008-01-02 08:31:28 UTC (rev 8975)
@@ -58,39 +58,43 @@
>>>
-
+> hc -help
hc - perform Hager & O'Connell flow computation
-This code can compute velocities, tractions, and geoid for simple
-density distributions and plate velocities using the semi-analytical
-approach of Hager & O'Connell (1981). This particular implementation
-illustrates one possible way to combine the HC solver routines.
-
+This code can compute velocities, tractions, and geoid for simple density distributions
+and plate velocities using the semi-analytical approach of Hager & O'Connell (1981).
+This particular implementation illustrates one possible way to combine the HC solver routines.
Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.
+This version by Thorsten Becker and Craig O'Neill
density anomaly options:
-dens name use name as a SH density anomaly model (dens.sh.dat)
- This expects density anomalies in % of PREM in DT format.
--ds additional density anomaly scaling factor (0.2)
+ All density anomalies are in units of 1% of PREM, all SH coefficients in Dahlen & Tromp convention.
+-dshs use the short, Becker & Boschi (2002) format for the SH density model (OFF)
+-ds density scaling factor (0.2)
+-dsf file read depth dependent density scaling from file
+ (overrides -ds, OFF), use pdens.py to edit
Earth model options:
--prem name set Earth model to name (/home/becker/progs/src/hc-svn/prem/prem.dat)
--vf name set viscosity structure filename to name (visc.dat)
+-prem name set Earth model to name (/home/walter/becker/progs/src/hc-svn/prem/prem.dat)
+-vf name viscosity structure filename (visc.dat), use pvisc.py to edit
This file is in non_dim_radius viscosity[Pas] format
boundary condition options:
--fs perform free slip computation, else no slip or plates (OFF)
--ns perform no slip computation, else try to read plate velocities (OFF)
--pvel name set surface velocity file to name (pvel.sh.dat)
- This file is based on a DT expansion of cm/yr velocity fields.
+-fs perform free slip computation (ON)
+-ns perform no slip computation (OFF)
+-pvel name set prescribed surface velocities from file name (OFF)
+ The file (e.g. pvel.sh.dat) is based on a DT expansion of cm/yr velocity fields.
solution procedure and I/O options:
-ng do not compute and print the geoid (OFF)
-pptsol print pol[6] and tor[2] solution vectors (OFF)
-px print the spatial solution to file (OFF)
--str compute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)
+-rtrac compute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)
+-htrac compute stt,stp,spp tractions [MPa] instead of velocities [cm/yr] (default: vel)
-v -vv -vvv: verbosity levels (0)
+
<<<
SPHERICAL HARMONICS FORMAT
@@ -159,10 +163,19 @@
Unformatted list of radii (radius of layer/Earth radius) and viscosity
(in Pas) values, reads until end of file. Values determine each layer
-viscosity upward until the next entry.
+viscosity upward until the next entry. Use the graphical tool
+pvisc.py to edit such files.
+2) Depth dependent density scaling file
+r_i d_i
+Format as for the viscosity file, but d_i are the depth-dependent
+scaling factors (this overridings -ds). Use the graphical tool
+pdens.py to edit such files.
+
+
+
OUTPUT FILES
After a regular run, file sol.bin will have the velocity solution
@@ -220,10 +233,15 @@
cat etopo5.ab | sh_syn > etopo5.127.dat
+SEATREE
+HC is a module of the Solid Earth Teaching and Research Environment
+(SEATREE) which provides a graphical user interface to flow
+computations and plotting.
+
+https://geosys.usc.edu/projects/seatree/wiki/
+
-
-
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/hc.h 2008-01-02 08:31:28 UTC (rev 8975)
@@ -125,10 +125,10 @@
hc_boolean compressible; /* compressibility following Panasyuk
& Steinberger */
hc_boolean free_slip; /* surface mechanical boundary condition */
- hc_boolean vel_bc_zero; /*
- if false, plate velocities, else no
- slip if free_slip is false as well
+ hc_boolean no_slip; /*
+ zero surface velocities
*/
+ hc_boolean platebc; /* plate velocity BC */
HC_PREC dens_anom_scale; /* default density anomaly scaling to
go from PREM percent traveltime
anomalies to density anomalies */
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/hc_init.c 2008-01-02 08:31:28 UTC (rev 8975)
@@ -20,11 +20,11 @@
*/
p->compressible = FALSE; /* compressibility following Panasyuk
& Steinberger */
- p->free_slip = FALSE; /* surface mechanical boundary condition */
- p->vel_bc_zero = FALSE; /*
- if false, plate velocities, else no
- slip if free_slip is false as well
- */
+ /* surface mechanical boundary condition */
+ p->free_slip = TRUE; /* free slip? */
+ p->no_slip = FALSE; /* no slip boundary condition? */
+ p->platebc = FALSE; /* plate velocities? */
+
p->compute_geoid = TRUE; /* compute the geoid? */
p->dens_anom_scale = 0.2; /* default density anomaly scaling to
go from PREM percent traveltime
@@ -46,7 +46,6 @@
*/
strncpy(p->visc_filename,HC_VISC_FILE,HC_CHAR_LENGTH);
strncpy(p->dens_filename,HC_DENS_SH_FILE,HC_CHAR_LENGTH);
- strncpy(p->pvel_filename,HC_PVEL_FILE,HC_CHAR_LENGTH);
strncpy(p->prem_model_filename,PREM_MODEL_FILE,HC_CHAR_LENGTH);
}
@@ -104,7 +103,11 @@
{
int dummy=0;
/* mechanical boundary condition */
- hc->free_slip = p->free_slip;
+ if(p->free_slip){
+ if(p->no_slip || p->platebc)
+ HC_ERROR("hc_init_main","free slip and no slip does not make sense");
+ hc->free_slip = TRUE;
+ }
/*
set the default expansion type, input expansions will be
converted
@@ -118,39 +121,66 @@
/*
initialize viscosity structure from file
*/
- hc_assign_viscosity(hc,HC_INIT_FROM_FILE,
- p->visc_filename,p->verbose);
- if(p->vel_bc_zero){
- /* no slip (zero velocity) surface boundary conditions */
+ hc_assign_viscosity(hc,HC_INIT_FROM_FILE,p->visc_filename,p->verbose);
+
+ if(p->no_slip && (!p->platebc)){
+ /*
+
+ no slip (zero velocity) surface boundary conditions
+
+ */
if(p->free_slip)
- HC_ERROR("hc_init","vel_bc_zero and free_slip doesn't make sense");
- /* read in the densities */
+ HC_ERROR("hc_init","no slip and free_slip doesn't make sense");
+ if(p->verbose)
+ fprintf(stderr,"hc_init: initializing for no slip\n");
+
+ /*
+ read in the densities first to determine L from the density expansion
+ */
hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
p->dens_filename,-1,FALSE,FALSE,p->verbose,p->read_short_dens_sh,
p->read_dens_scale_from_file,p->dens_scaling_filename);
- /* assign all zeroes up to the lmax of the density expansion */
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,p->vel_bc_zero,hc->dens_anom[0].lmax,
+ /*
+ assign all zeroes up to the lmax of the density expansion
+ */
+ hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,TRUE,hc->dens_anom[0].lmax,
FALSE,p->verbose);
+ }else if(p->platebc){
+ /*
+
+ surface velocities
+
+ */
+ if(p->free_slip)
+ HC_ERROR("hc_init","plate boundary condition and free_slip doesn't make sense");
+ if(p->verbose)
+ fprintf(stderr,"hc_init: initializing for surface velocities\n");
+ /*
+ read in velocities, which will determine the solution lmax
+ */
+ hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,FALSE,dummy,FALSE,p->verbose);
+ /* then read in the density anomalies */
+ hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,hc->pvel[0].lmax,
+ FALSE,FALSE,p->verbose,p->read_short_dens_sh,
+ p->read_dens_scale_from_file,p->dens_scaling_filename);
+ }else if(p->free_slip){
+ /*
+
+ free slip
+
+ */
+ if(p->no_slip)
+ HC_ERROR("hc_init","no slip and free slip does not make sense");
+ if(p->verbose)
+ fprintf(stderr,"hc_init: initializing for free-slip\n");
+ /* read in the density fields */
+ hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,-1,FALSE,FALSE,
+ p->verbose,p->read_short_dens_sh,
+ p->read_dens_scale_from_file,p->dens_scaling_filename);
}else{
- /* presribed surface velocities */
- if(!p->free_slip){
- /* read in velocities, which will determine the solution lmax */
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,
- p->pvel_filename,
- p->vel_bc_zero,
- dummy,FALSE,p->verbose);
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,hc->pvel[0].lmax,
- FALSE,FALSE,p->verbose,p->read_short_dens_sh,
- p->read_dens_scale_from_file,p->dens_scaling_filename);
- }else{
- if(p->verbose)
- fprintf(stderr,"hc_init: initializing for free-slip\n");
- /* read in the density fields */
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,-1,FALSE,FALSE,
- p->verbose,p->read_short_dens_sh,
- p->read_dens_scale_from_file,p->dens_scaling_filename);
- }
+ HC_ERROR("hc_init","boundary condition logic error");
}
+
/*
phase boundaries, if any
*/
@@ -274,25 +304,24 @@
fprintf(stderr,"-dshs\t\tuse the short, Becker & Boschi (2002) format for the SH density model (%s)\n",
hc_name_boolean(p->read_short_dens_sh));
- fprintf(stderr,"-ds\t\tadditional density anomaly scaling factor (%g)\n",
+ fprintf(stderr,"-ds\t\tdensity scaling factor (%g)\n",
p->dens_anom_scale);
- fprintf(stderr,"-dsf\tfile\tread depth dependent density scaling from file (overrides -ds, %s)\n\n",
+ fprintf(stderr,"-dsf\tfile\tread depth dependent density scaling from file\n");
+ fprintf(stderr,"\t\t(overrides -ds, %s), use pdens.py to edit\n\n",
hc_name_boolean(p->read_dens_scale_from_file));
fprintf(stderr,"Earth model options:\n");
fprintf(stderr,"-prem\tname\tset Earth model to name (%s)\n",
p->prem_model_filename);
- fprintf(stderr,"-vf\tname\tset viscosity structure filename to name (%s)\n",
+ fprintf(stderr,"-vf\tname\tviscosity structure filename (%s), use pvisc.py to edit\n",
p->visc_filename);
fprintf(stderr,"\t\tThis file is in non_dim_radius viscosity[Pas] format\n\n");
fprintf(stderr,"boundary condition options:\n");
- fprintf(stderr,"-fs\t\tperform free slip computation, else no slip or plates (%s)\n",
- hc_name_boolean(p->free_slip));
- fprintf(stderr,"-ns\t\tperform no slip computation, else try to read plate velocities (%s)\n",
- hc_name_boolean(p->vel_bc_zero));
- fprintf(stderr,"-pvel\tname\tset surface velocity file to name (%s)\n",
- p->pvel_filename);
- fprintf(stderr,"\t\tThis file is based on a DT expansion of cm/yr velocity fields.\n\n");
+ fprintf(stderr,"-fs\t\tperform free slip computation (%s)\n",hc_name_boolean(p->free_slip));
+ fprintf(stderr,"-ns\t\tperform no slip computation (%s)\n",hc_name_boolean(p->no_slip));
+ fprintf(stderr,"-pvel\tname\tset prescribed surface velocities from file name (%s)\n",
+ hc_name_boolean(p->platebc));
+ fprintf(stderr,"\t\tThe file (e.g. %s) is based on a DT expansion of cm/yr velocity fields.\n\n",HC_PVEL_FILE);
fprintf(stderr,"solution procedure and I/O options:\n");
fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%s)\n",
@@ -311,9 +340,9 @@
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],HC_FLT_FORMAT,&p->dens_anom_scale);
}else if(strcmp(argv[i],"-fs")==0){ /* free slip flag */
- hc_toggle_boolean(&p->free_slip);
+ p->free_slip = TRUE;p->no_slip = FALSE;
}else if(strcmp(argv[i],"-ns")==0){ /* no slip flag */
- hc_toggle_boolean(&p->vel_bc_zero);
+ p->no_slip = TRUE;p->free_slip = FALSE;
}else if(strcmp(argv[i],"-px")==0){ /* print spatial solution? */
hc_toggle_boolean(&p->print_spatial);
}else if(strcmp(argv[i],"-dshs")==0){ /* use short format? */
@@ -331,9 +360,10 @@
}else if(strcmp(argv[i],"-prem")==0){ /* PREM filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->prem_model_filename,argv[i],HC_CHAR_LENGTH);
- }else if(strcmp(argv[i],"-pvel")==0){ /* velocity filename */
+ }else if(strcmp(argv[i],"-pvel")==0){ /* velocity filename, this will switch off free slip */
hc_advance_argument(&i,argc,argv);
strncpy(p->pvel_filename,argv[i],HC_CHAR_LENGTH);
+ p->platebc = TRUE;p->no_slip = TRUE;p->free_slip = FALSE;
}else if(strcmp(argv[i],"-dens")==0){ /* density filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_filename,argv[i],HC_CHAR_LENGTH);
@@ -782,9 +812,8 @@
int i;
if(depth_dependent){
i=0;
- while((i<n) && (rd[i] < r)){
+ while((i<n) && (rd[i] < r))
i++;
- }
i--;
return sd[i];
}else{
@@ -829,11 +858,8 @@
*/
-void hc_assign_plate_velocities(struct hcs *hc,int mode,
- char *filename,
- hc_boolean vel_bc_zero,int lmax,
- hc_boolean pvel_in_binary,
- hc_boolean verbose)
+void hc_assign_plate_velocities(struct hcs *hc,int mode, char *filename,hc_boolean vel_bc_zero,
+ int lmax,hc_boolean pvel_in_binary,hc_boolean verbose)
{
int type,shps,ilayer,nset,ivec;
HC_PREC zlabel,vfac[2],t10[2],t11[2];
Added: mc/1D/hc/trunk/manipulate_hc.py
===================================================================
--- mc/1D/hc/trunk/manipulate_hc.py 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/manipulate_hc.py 2008-01-02 08:31:28 UTC (rev 8975)
@@ -0,0 +1,382 @@
+#
+# manipulate a HC type data profile
+#
+import math
+import pylab as p
+import matplotlib
+from matplotlib.backends.backend_gtk import FigureCanvasGTK, NavigationToolbar
+matplotlib.use('GTK')
+
+
+class ManipulateXYData:
+ """
+
+ handles x y data that specifies HC-type profiles
+
+ use like:
+
+ >>>
+
+ mp = ManipulateXYData(filename,mode)
+ p.connect('button_press_event', mp.on_click)
+ p.connect('button_release_event', mp.on_release)
+
+ <<<
+
+ INPUT
+
+ filename: data name
+
+ mode: 1: data are viscosity
+ 2: data are density scaling factors
+
+
+ xtol amd ytol are relative tolerances
+
+ inspired by http://www.scipy.org/Cookbook/Matplotlib/Interactive_Plotting
+
+ """
+
+ # initialize class
+ def __init__(self, filename, mode, xtol = None, ytol = None):
+ if xtol == None:
+ xtol = 0.1
+ if ytol == None:
+ ytol = 0.1
+ self.xtol = xtol
+ self.ytol = ytol
+
+
+ self.visc_norm = 1.e21 # some constants
+ self.radius_km = 6371.
+ self.cmb_km = 2891.
+
+
+ self.figure = p.figure()
+ self.axis = self.figure.add_subplot(111)
+
+ #
+ # 1: viscosity
+ # 2: density
+ self.plot_mode = mode
+ #
+ # read data
+ self.read_data(filename,mode)
+ #
+ # convert to plotting format
+ self.convert_data(self.plot_mode,False)
+
+ self.datax0 = self.datax # copy for restore
+ self.datay0 = self.datay
+
+ self.moving = -1 # if point are being move
+
+
+ self.zlabels = [300,660,1750] # for plot
+ #
+ self.verbose = True # progress output
+
+ self.xl = []
+ self.yl = []
+
+ # start a plot
+ self.redraw_plot()
+
+
+ def distance(self, x1, x2, y1, y2):
+ """
+ return the distance between two points
+ """
+ return(math.sqrt( (x1 - x2)**2 + (y1 - y2)**2 ))
+
+
+ def __call__(self, event):
+ #
+ # get the x and y data coords
+ #
+ x, y = event.xdata, event.ydata
+
+ if event.inaxes:
+ print 'generic call?: data coords', x,y
+
+ def on_click(self, event):
+ #
+ # mouse click event
+ #
+ if event.inaxes: # within region
+ # data coordinates
+ x, y = event.xdata, event.ydata
+ if self.axis == event.inaxes:
+ #
+ # look for close points
+ #
+ cps = []
+ i=0
+ data_cont = zip(self.datax,self.datay) # reformat
+ for xd,yd in data_cont: # compute distances for
+ # those points that are
+ # within range compute
+ # tolerance
+ if xd == 0: # compute tolerances
+ xts = 0.5
+ else:
+ xts = abs(xd)*self.xtol
+ if yd == 0:
+ yts = 0.5
+ else:
+ yts = abs(yd)*self.ytol
+ #
+ # if close, compute distance
+ if (abs(x-xd) < xts) and (abs(y-yd) < yts) :
+ cps.append( (self.distance(xd,x,yd,y),xd,yd,i) )
+ i=i+1
+ if cps: # if we found some point close enough, sort them and use the closest
+ cps.sort()
+ dist, xd, yd, i = cps[0] # closest
+
+ if event.button == 2: # center mouse click: add point to list
+ if not cps or dist > 1:
+ if self.verbose:
+ print 'adding data point %7.2f, %7.2f ' % (x, y)
+ self.datax.append(x)
+ self.datay.append(y)
+ self.redraw_plot()
+ else:
+ if self.verbose:
+ print 'there is already a point at %7.2f, %7.2f ' % (x, y)
+ else:
+ #
+ # left or right
+ #
+ if cps:
+ if event.button == 1:
+ # left mouse button, move this data point
+ self.moving = i
+ if self.verbose:
+ print 'moving data point %5i ' % i, 'from %7.2f, %7.2f ' % (xd, yd)
+ elif event.button == 3:
+ # right click: removing this data point
+ if self.verbose:
+ print 'removing data point %7.2f, %7.2f ' % (self.datax[i],self.datay[i])
+ del self.datax[i]
+ del self.datay[i]
+ self.redraw_plot()
+ else:
+ if self.verbose:
+ print 'did not find data close to click %7.2f, %7.2f' % (x,y)
+
+
+ def on_release(self, event):
+ # mouse has been released
+ if event.inaxes:
+ xd, yd = event.xdata, event.ydata
+ if self.axis == event.inaxes:
+ if self.moving > -1: # are we actually moving a point?
+ if self.verbose:
+ print 'assigning %7.2f, %7.2f to data point %5i' % (xd, yd, self.moving)
+ i=0;xn,yn=[],[]
+ data_cont=zip(self.datax,self.datay)
+ # this could be dealt with smarter
+ self.datax, self.datay = [], []
+ for x,y in data_cont: # replace the self.moving-th point with the current location
+ if i==self.moving:
+ self.datax.append(xd);self.datay.append(yd)
+ else:
+ self.datax.append(x);self.datay.append(y)
+ i+=1
+ self.redraw_plot()
+ self.moving = -1 # reset
+
+ def redraw_plot(self): # refresh the plot
+ """
+
+ redraw a plot
+
+ """
+ self.sortlevels() # sort data and get layer entries
+
+ # get the figure handle
+ p.axes(self.axis)
+ self.axis.clear()
+ if self.plot_mode == 1: # viscosity
+ xmin,xmax = 1e-3,1e3
+ if min(self.datax) < xmin:
+ xmin /= 10.
+ if max(self.datax) > xmax:
+ xmax *= 10.
+ self.axis.semilogx(self.datax,self.datay,'o') # plot actual profile
+ self.axis.semilogx(self.xl,self.yl,linewidth=3,color='red') # plot layers
+ self.axis.set_xlabel('viscosity [1e21 Pas]')
+ self.add_pmantle_ornaments(xmin,xmax,True)
+
+ elif self.plot_mode == 2: # density scaling factor
+ xmin,xmax = -0.1,0.4
+ if min(self.datax) < xmin:
+ xmin = self.datax *0.8
+ if max(self.datax) > xmax:
+ xmax = self.datay *1.2
+ self.axis.plot(self.datax,self.datay,'o') # plot actual profile
+ self.axis.plot(self.xl,self.yl,linewidth=3,color='blue')
+ self.axis.set_title('left mouse: move center: add right: remove point')
+ self.axis.set_xlabel('scale factor')
+ self.add_pmantle_ornaments(xmin,xmax,False)
+# what is the renderer?
+# self.axis.draw('GTKAgg')
+ p.draw()
+
+
+ def add_pmantle_ornaments(self,xmin,xmax,uselogx):
+ """
+ add ornaments typical for the earth's mantle to the plot
+ """
+ self.axis.grid(True)
+ self.axis.set_title('left mouse: move center: add right: remove point')
+ self.axis.set_ylabel('depth [km]')
+ x = [xmin,xmax];
+
+ if self.plot_mode == 1:
+ xoff = xmin*2.5
+ else:
+ xoff = 0.025*(xmax-xmin)
+
+ for z in self.zlabels: # add a few labels
+ y=[-z,-z];
+ self.axis.text(xmin+xoff,-z+10.,str(z)+' km',fontstyle='italic')
+ if uselogx:
+ self.axis.semilogx(x,y,linewidth=2,color='black',linestyle='--')
+ else:
+ self.axis.plot(x,y,linewidth=2,color='black',linestyle='--')
+ self.axis.set_ylim([-self.cmb_km, 0])
+ self.axis.set_xlim([xmin,xmax])
+
+
+
+ def reset_data(self,event):
+ if self.verbose:
+ print 'resetting to original data'
+ self.datax = self.datax0
+ self.datay = self.datay0
+ self.redraw_plot()
+
+ def sortlevels(self):
+ """
+
+ sort through a list of weirdly formatted viscosity file
+ values and add data point to make a plot look nice
+
+ also, assign layer plot data
+
+ """
+ # sort the z and eta vectors
+ data = []
+ for zl, el in zip(self.datay, self.datax):
+ data.append((zl, el))
+ data.sort();
+ z,eta = [], []
+ for zl,el in data:
+ z.append(zl); eta.append(el)
+
+ zn, en =[], []
+ n = len(z)
+ if n:
+ if z[0] > -self.cmb_km:
+ zn.append(-self.cmb_km)
+ en.append(eta[0])
+ i=0
+ while i < n:
+ zn.append(z[i])
+ en.append(eta[i])
+ i += 1
+ if n and z[n-1] < 0:
+ zn.append(0)
+ en.append(eta[n-1])
+ self.datax = en
+ self.datay = zn
+
+
+ # convert the point-based data to one that can be plotted as
+ # layers
+ data = []
+ i=0; n = len(self.datax);
+ while i < n:
+ if i > 0 and self.datay[i] != self.datay[i-1]:
+ data.append((self.datay[i],self.datax[i-1]))
+ data.append((self.datay[i],self.datax[i]))
+ i += 1
+ self.xl,self.yl = [],[]
+ for yl,xl in data:
+ self.xl.append(xl)
+ self.yl.append(yl)
+
+ def read_data(self,filename,mode):
+ """
+
+ read HC profile data from file and return datax, datay
+
+ """
+ self.datax,self.datay = [],[]
+ f=open(filename,'r')
+ for line in f:
+ val = line.split()
+ if len(val) != 2:
+ print 'error file ', filename, ' appears to be in wrong format'
+ print 'expecting'
+ if self.mode == 1:
+ print 'radius[non_dim] viscosity[Pas]'
+ elif self.mode == 2:
+ print 'radius[non_dim] density_scale'
+ else:
+ print 'unknown'
+ exit()
+ self.datax.append(val[0])
+ self.datay.append(val[1])
+ f.close()
+
+
+
+ def convert_data(self,mode,reverse):
+ """ convert input data to plotting format and vice versa """
+ tmpx, tmpy = [],[]
+ i=0
+ for i in range(0,len(self.datax)):
+ if mode == 1: # viscosity
+ if not reverse:
+ tmpx.append(float(self.datay[i])/self.visc_norm)
+ tmpy.append(-(1.-float(self.datax[i]))*self.radius_km)
+ else:
+ tmpy.append(float(self.datax[i])*self.visc_norm)
+ tmpx.append((self.radius_km+float(self.datay[i]))/self.radius_km)
+ elif mode == 2: # density
+ if not reverse:
+ tmpx.append(float(self.datay[i]))
+ tmpy.append(-(1.-float(self.datax[i]))*self.radius_km)
+ else:
+ tmpy.append(float(self.datax[i]))
+ tmpx.append((self.radius_km+float(self.datay[i]))/self.radius_km)
+ self.datax,self.datay = tmpx,tmpy
+
+
+ def save_and_exit(self,event):
+ if self.verbose:
+ print 'saving modified data'
+ filename = 'tmp.dat'
+ if self.plot_mode == 1:
+ print 'saving modified viscosity profile data to ', filename
+ elif self.plot_mode == 2:
+ print 'saving modified density profile data to ', filename
+
+ #
+ # convert data back
+ self.convert_data(self.plot_mode,True)
+ f=open(filename,'w')
+ i=0
+ for i in range(0,len(self.datax)):
+ ostring = "%8.5f\t%12.7e\n" % (self.datax[i], self.datay[i])
+ f.writelines(ostring)
+ f.close()
+ p.close(self.figure)
+
+ def exit(self,event):
+ if self.verbose:
+ print 'exiting without saving'
+ p.close(self.figure)
Added: mc/1D/hc/trunk/pdens.py
===================================================================
--- mc/1D/hc/trunk/pdens.py 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/pdens.py 2008-01-02 08:31:28 UTC (rev 8975)
@@ -0,0 +1,39 @@
+#! /usr/bin/env python
+#
+# plot and mofiy a HC density scaling file
+#
+from optparse import OptionParser
+import math
+import pylab as p
+import matplotlib
+from matplotlib.backends.backend_gtk import FigureCanvasGTK, NavigationToolbar
+matplotlib.use('GTK')
+from manipulate_hc import *
+
+
+parser = OptionParser()
+parser.add_option("--df", "--density-scaling-file", type="string", dest="filename", default="dscale.dat", \
+ help="HC type density scaling file to display and modify")
+(options, args) = parser.parse_args()
+filename = options.filename
+print 'using ', filename, ' as HC density scaling profile '
+
+mp = ManipulateXYData(filename,2)
+#
+# connect main data window plot handling routines with mouse click procedures
+p.connect('button_press_event', mp.on_click)
+p.connect('button_release_event', mp.on_release)
+# add some GUI neutral buttons
+reset_button = matplotlib.widgets.Button(p.axes([0.85, 0.025, 0.1, 0.04]), 'Reset', hovercolor='gray')
+cancel_button = matplotlib.widgets.Button(p.axes([0.75, 0.025, 0.1, 0.04]), 'Cancel', hovercolor='gray')
+done_button = matplotlib.widgets.Button(p.axes([0.65, 0.025, 0.1, 0.04]), 'Save', hovercolor='gray')
+# button handling routines
+reset_button.on_clicked(mp.reset_data)
+done_button.on_clicked(mp.save_and_exit)
+cancel_button.on_clicked(mp.exit)
+
+p.show()
+
+
+
+
Property changes on: mc/1D/hc/trunk/pdens.py
___________________________________________________________________
Name: svn:executable
+ *
Added: mc/1D/hc/trunk/pvisc.py
===================================================================
--- mc/1D/hc/trunk/pvisc.py 2007-12-31 02:25:00 UTC (rev 8974)
+++ mc/1D/hc/trunk/pvisc.py 2008-01-02 08:31:28 UTC (rev 8975)
@@ -0,0 +1,41 @@
+#! /usr/bin/env python
+#
+# plot and mofiy a HC viscosity file
+#
+from optparse import OptionParser
+import math
+import pylab as p
+import matplotlib
+from matplotlib.backends.backend_gtk import FigureCanvasGTK, NavigationToolbar
+matplotlib.use('GTK')
+from manipulate_hc import *
+
+
+parser = OptionParser()
+parser.add_option("--vf", "--viscosity-file", type="string", dest="filename", default="visc.C", \
+ help="HC type viscosity file to display and modify")
+(options, args) = parser.parse_args()
+filename = options.filename
+
+print 'using ', filename, ' as HC viscosity profile '
+
+
+mp = ManipulateXYData(filename,1)
+#
+# connect main data window plot handling routines with mouse click procedures
+p.connect('button_press_event', mp.on_click)
+p.connect('button_release_event', mp.on_release)
+# add some GUI neutral buttons
+reset_button = matplotlib.widgets.Button(p.axes([0.85, 0.025, 0.1, 0.04]), 'Reset', hovercolor='gray')
+cancel_button = matplotlib.widgets.Button(p.axes([0.75, 0.025, 0.1, 0.04]), 'Cancel', hovercolor='gray')
+done_button = matplotlib.widgets.Button(p.axes([0.65, 0.025, 0.1, 0.04]), 'Save', hovercolor='gray')
+# button handling routines
+reset_button.on_clicked(mp.reset_data)
+done_button.on_clicked(mp.save_and_exit)
+cancel_button.on_clicked(mp.exit)
+
+p.show()
+
+
+
+
Property changes on: mc/1D/hc/trunk/pvisc.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the cig-commits
mailing list