|
Matlab/Octave tools for geophysical studiesFrançois Beauducel |
| Model | Description | Parameters | Outputs | References |
| |||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Space | Source | Medium | ||||||||||
| Mogi | Point source in elastic half-space Approximation for sphere of radius a << f |
r | f, ΔV f, ΔP, a |
ν E, ν or μ |
displacements, tilt, strains at surface | [Anderson, 1936] [Mogi, 1958] |
= 1.5 | |||||
| Sun | Penny-shaped crack in elastic half-space Approximation for h/a >> 1 |
r | h, a, ΔV h, a, ΔP |
E, ν |
displacements at surface | [Sun, 1969] | ||||||
| Okada | Rectangular dislocation in elastic half-space | x,y | d, δ, θ, W, L R, USLIP, UOPEN |
ν | displacements, tilt, strains at surface | [Okada, 1985] |
|
= 1 (sill) = 0.75 (dyke) |
||||
| Okubo | Rectangular dislocation in elastic half-space | x,y | d, δ, θ, W, L R, USLIP, UOPEN |
ν, ρ, ρ' | gravity and elevation change at surface | [Okubo, 1992] | ||||||
| r | Radial distance from source |
|---|---|
| x | X-coordinate from source |
| y | Y-coordinate from source |
| f, h | Focal depth of source |
| a | Radius of source |
| d | Depth of top-edge fault |
| W | Width of fault |
| L | Length of fault |
| δ | Dip-angle of fault |
| θ | Strike-angle of fault (azimut) |
| ΔV | Volume variation of source |
| ΔP | Pressure variation of source |
| R | Rake-angle of dislocation |
| USLIP | Slip-dislocation (in Rake-angle direction) |
| UOPEN | Tensile-dislocation (opening) |
| μ, G | Rigidity of medium (Lamé's second parameter) |
| λ | Lamé's first parameter |
| E | Elasticity of medium (Young's modulus) |
| ν | Poisson's ratio of medium |
| ρ | Density of medium |
| ρ' | Density of cavity-filling matter |
Some conversion formulas between elastic parameters λ, μ, E and ν:
| ||||||||||||
| ||||||||||||
| ||||||||||||
| ||||||||||||
|
 |
The Mogi [1958] model calculates analytical solution for surface deformation due to a point source in an elastic half-space. This model is widely used to simulate ground deformation produced by local perturbation like volcanic magma chamber. It computes displacements, tilt and strain in a polar space, due to a volume variation at depth [Anderson, 1936] or an isotropic pressure variation in a spherical source [Mogi, 1958]. The proposed Matlab script is a literal transcription of Mogi's simple equations, extended to non-isotropic medium (Poisson's ratio different from 0.25). All parameters can be vectorized. See help for syntax, and script comments for details. MOGI Mogi's model (point source in elastic half-space). [Ur,Uz,τ,σr,σt] = MOGI(R,F,ΔV,ν) or MOGI(R,F,A,ΔP,E,ν) computes radial and vertical displacements Ur and Uz, ground tilt τ, radial and tangential strain σr and σt on surface, at a radial distance R from the top of the source due to a hydrostatic pressure inside a sphere of radius A at depth F, in a homogeneous, semi-infinite elastic body and approximation for A << (center of dilatation). Formula by Anderson [1936] and Mogi [1958]. MOGI(R,F,ΔV) and MOGI(R,F,A,μ,ΔP) are also allowed for compatibility (Mogi's original equation considers an isotropic material with Lamé constants equal, i.e., λ = μ, ν = 0.25). Input variables are: F: depth of the center of the sphere from the surface, ΔV: volumetric change of the sphere, A: radius of the sphere, ΔP: hydrostatic pressure change in the sphere, E: elasticity (Young's modulus), ν: Poisson's ratio, μ: rigidity (Lamé constant in case of isotropic material). Notes: - Equations are all vectorized, so variables R,F,V,A,μ and P can be vectors or N-D matrix of the same size and any of them can be scalars; then outputs will be also vector or matrix. - Convention: Uz > 0 = UP, F is depth so in -Z direction. - Units should be constistent, e.g.: R, F, A, Ur and Uz in m imply ΔV in m³; E, μ and P in Pa; τ in rad, σr, σt and ν dimensionless. Example for a 3-D plot of exagerated deformed surface due to a 1-bar overpressure in a 10-cm radius sphere at 1-m depth in rock: [x,y] = meshgrid(-3:.1:3); [th,rho] = cart2pol(x,y); [ur,uz] = mogi(rho,1,0.1,1e5,10e9,0.25); [ux,uy] = pol2cart(th,ur); ps = 1e8; surf(x+ux*ps,y+uy*ps,uz*ps), axis equal, light Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 1997 Updated: 2012-04-05 References: Anderson, E.M., Dynamics of the formation of cone-sheets, ring-dikes, and cauldron-subsidences, Proc. R. Soc. Edinburgh, 56, 128-157, 1936. Mogi, K., Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them, Bull. Earthquake Res. Inst. Univ. Tokyo, 36, 99-134, 1958. Acknowledgments: Raphaël Grandin, Benoît Taisne Download mogi.m at Matlab Central file exchange (Copyright 1997-2012 F. Beauducel / BSD License). |
 |
The Sun [1969] model calculates analytical solution for surface deformation due to hydrostatic pressure inside a horizontal circular fracture (“penny-shaped”) in an elastic half-space. The proposed Matlab script is a literal transcription of the Sun's paper equations. The equations are also vectorized for all input parameters. See help for syntax, and script comments for details. SUN69 Deformation from penny-shaped crack in elastic half-space. [Ur,Uz] = SUN69(R,H,A,ΔV) or SUN69(R,H,A,ΔP,E,ν) computes radial and vertical displacements Ur and Uz on the free surface, due to a horizontal circular fracture formed in a semi-infinite elastic medium, with following variables: R: radial distance of observation, H: depth of the center of the source from the surface, A: radius of the source with the hydrostatic pressure, ΔV: volume of injected material, ΔP: change of the hydrostatic pressure in the crack. E: Young's modulus, ν: Poisson's ratio (default is 0.25 for isotropic medium). [Ur,Uz,B] = SUN69(...) returns also the maximum separation of fracture (vertical displacement) B. Equations from Sun [1969], with approximation H/A >> 1. If H/A > 2, error is about 2 to 3%; if H/A > 5, solution is almost perfect. Notes: - Equations are all vectorized, so variables R,H,A,ΔV or ΔP can be scalar or any of them vector or matrix, then outputs will have the same size. - Convention: Uz > 0 = UP, H is depth so in -Z direction. - Units should be constistent, e.g.: R, H, A, Ur and Uz in m imply V in m³; optional E and ΔP in Pa, ν dimensionless. Example for a 3D plot of exagerated deformed surface: [x,y] = meshgrid(-3:.1:3); [th,rho] = cart2pol(x,y); [ur,uz] = sun69(rho,1,.5,1e6,10e9,0.25); [ux,uy] = pol2cart(th,ur); ps = 5e4; surf(x+ux*ps,y+uy*ps,uz*ps), axis equal, light Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris, 2009. References: Sun, R. J. (1969). Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium, J. Geophys. Res., 74, 5995-6011. Download sun69.m at Matlab Central file exchange (Copyright 2010 F. Beauducel / BSD License). |
 |
The Okada [1985] model calculates analytical solution for surface deformation due to shear and tensile faults in an elastic half-space. This model is widely used to simulate ground deformation produced by local perturbation like tectonic faults (earthquakes) or volcanic dykes (magmatic intrusion). Given rectangular fault geometry (length, width, depth, strike, dip) and 3-component dislocation amplitude (rake, slip and open), it computes the displacements, tilts and strains at the free-surface. The proposed Matlab script is a literal transcription of the Okada's equations, except that it is transposed in a geographical referential (East, North, Up), where the fault is defined by a strike angle relative to the North, and dislocation parameters are given by: rake, slip and opening (instead of U1, U2, U3), following Aki & Richards [1980] definition. All coordinates and depth are relative to fault centroid. Lamé's constants λ and μ are replaced by Poisson's ratio ν (with a default value of 0.25 for isotropic medium), since the equations are independent of other elastic parameters. The equations are also vectorized for (x,y) coordinates and all input parameters except dip angle. To check the consistency of numerical calculations, run the script okada85_checklist.m, a transcription of table 2 cases 2, 3, and 4 checklist from [Okada, 1985] paper. See help for further details, syntax, example, and script comments for technical details. OKADA85 Surface deformation due to a finite rectangular source. [uE,uN,uZ,uZE,uZN,uNN,uNE,uEN,uEE] = OKADA85(... E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN) computes displacements, tilts and strains at the surface of an elastic half-space, due to a dislocation defined by RAKE, SLIP, and OPEN on a rectangular fault defined by orientation STRIKE and DIP, and size LENGTH and WIDTH. The fault centroid is located (0,0,-DEPTH). E,N : coordinates of observation points in a geographic referential (East,North,Up) relative to fault centroid (units are described below) DEPTH : depth of the fault centroid (DEPTH > 0) STRIKE : fault trace direction (0 to 360° relative to North), defined so that the fault dips to the right side of the trace DIP : angle between the fault and a horizontal plane (0 to 90°) LENGTH : fault length in the STRIKE direction (LENGTH > 0) WIDTH : fault width in the DIP direction (WIDTH > 0) RAKE : direction the hanging wall moves during rupture, measured relative to the fault STRIKE (-180 to 180°) SLIP : dislocation in RAKE direction (length unit) OPEN : dislocation in tensile component (same unit as SLIP) returns the following variables (same matrix size as E and N): uN,uE,uZ : displacements (unit of SLIP and OPEN) uZE,uZN : tilts (in rad * FACTOR) uNN,uNE,uEN,uEE : horizontal strains POSITIVE = COMPRESSION (unit of FACTOR) Length unit consistency: E, N, DEPTH, LENGTH, and WIDTH must have the same unit (e.g. km) which can be different from that of SLIP and OPEN (e.g. m) but with a possible FACTOR on tilt and strain results (in this case, an amplification of km/m = 1000). To have FACTOR = 1 (tilt in radians and correct strain unit), use the same length unit for all aforesaid variables. [...] = OKADA85(...,NU) specifies Poisson's ratio NU (default is 0.25 for an isotropic medium). Formulas and notations from Okada [1985] solution excepted for strain convention (here positive strain means compression), and for the fault parameters after Aki & Richards [1980], e.g.: DIP=90, RAKE=0 : left lateral (senestral) strike slip DIP=90, RAKE=180 : right lateral (dextral) strike slip DIP=70, RAKE=90 : reverse fault DIP=70, RAKE=-90 : normal fault It is also possible to produce partial outputs, with following syntax: [uE,uN,uZ] = OKADA85(...) for displacements only; [uE,uN,uZ,uZE,uZN] = OKADA85(...) for displacements and tilts; [uE,uN,uZ,uNN,uNE,uEN,uEE] = OKADA85(...) for displacements and strains; [uZE,uZN] = OKADA85(...) for tilts only; [uZE,uZN,uNN,uNE,uEN,uEE] = OKADA85(...) for tilts and strains; [uNN,uNE,uEN,uEE] = OKADA85(...) for strains only. Note that vertical strain components can be obtained with following equations: uNZ = -uZN; uEZ = -uZE; uZZ = -(uEE + uNN)*NU/(1-NU); [...] = OKADA85(...,'plot') or OKADA85(...) without output argument produces a 3-D figure with fault geometry and dislocation at scale (if all of the fault parameters are scalar). Equations are all vectorized excepted for argument DIP which must be a scalar (beacause of a singularity in Okada's equations); all other arguments can be scalar or matrix of the same size. Example: [E,N] = meshgrid(linspace(-10,10,50)); [uE,uN,uZ] = okada85(E,N,2,30,70,5,3,-45,1,1,'plot'); figure, surf(E,N,uN) considers a 5x3 fault at depth 2, N30° strike, 70° dip, and unit dislocation in all directions (reverse, senestral and open). Displacements are computed on a regular grid from -10 to 10, and North displacements are plotted as a surface. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 1997 Updated: 2014-05-24 References: Aki K., and P. G. Richards, Quantitative seismology, Freemann & Co, New York, 1980. Okada Y., Surface deformation due to shear and tensile faults in a half-space, Bull. Seismol. Soc. Am., 75:4, 1135-1154, 1985. Acknowledgments: Dmitry Nicolsky, Qian Yao, Halldor Geirsson. : new version on May 24, 2014 with bug corrections. Update your file! Download okada85.m at Matlab Central file exchange (Copyright 1997-2014 F. Beauducel / BSD License).
|
 |
The Okubo [1992] model calculates analytical solution for gravity changes due to shear and tensile faults in an elastic half-space (as Okada [1985] does for deformations). This model can be used to simulate gravity anomalies produced by local perturbation like tectonic faults (earthquakes) or volcanic dykes (magmatic intrusion). Given rectangular fault geometry (length, width, depth, strike, dip) and 3-component dislocation amplitude (rake, slip and open), it computes the total gravity change (potential, cavity filling and free-air) at the free-surface. The proposed Matlab script is a literal transcription of the Okubo's equations, except for the fault geometry that it is transposed in a geographical referential (East, North, Up), where the fault is defined by a strike angle relative to the North, and dislocation parameters are given by: rake, slip and opening (instead of U1, U2, U3), following Aki & Richards [1980] definition. All coordinates and depth are relative to fault centroid. The equations are also vectorized for (x,y) coordinates and all parameters except dip angle. Type "doc okubo92" for help, syntax and example, and see script comments for details. OKUBO92 Surface gravity change due to a finite rectangular source. [dG,dH] = OKUBO92(E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN,RHO,RHOP) computes total gravity change at the free surface of an elastic half-space, due to a dislocation defined by RAKE, SLIP, and OPEN on a rectangular fault defined by orientation STRIKE and DIP, and size LENGTH and WIDTH. The fault centroid is located (0,0,-DEPTH). E,N : coordinates of observation points in a geographic referential (East,North,Up) relative to fault centroid (units are described below) DEPTH : depth of the fault centroid (DEPTH > 0) STRIKE : fault trace direction (0 to 360° relative to North), defined so that the fault dips to the right side of the trace DIP : angle between the fault and a horizontal plane (0 to 90°) LENGTH : fault length in the STRIKE direction (LENGTH > 0) WIDTH : fault width in the DIP direction (WIDTH > 0) RAKE : direction the hanging wall moves during rupture, measured relative to the fault STRIKE (-180 to 180°). SLIP : dislocation in RAKE direction (in m) OPEN : dislocation in tensile component (in m) RHO : density of the medium (in kg/m^3) RHOP : optional density of the cavity-filling matter (in kg/m^3); takes RHO value if not specified. returns the following variables (same matrix size as E and N): dG : total gravity change (in m/s^2) dH : elevation changes as given by Okada [1985] model (in m) Length unit consistency: E, N, DEPTH, LENGTH, and WIDTH must have the same unit (e.g. km) which can be different from that of SLIP and OPEN (e.g. m to to produce dG in m/s^2 and dH in m). [...] = OKUBO92(...,RHOP,BETA,NU) specifies gravity vertical gradient BETA (default is 0.309e-5 s^-2) and Poisson's ratio NU (default is 0.25 for an isotropic medium). Formulas and notations from Okubo [1992] solution excepted for the fault geometry parameters after Aki & Richards [1980], e.g.: DIP=90, RAKE=0 : left lateral (senestral) strike slip DIP=90, RAKE=180 : right lateral (dextral) strike slip DIP=70, RAKE=90 : reverse fault DIP=70, RAKE=-90 : normal fault Equations are all vectorized excepted for argument DIP which must be a scalar; all other arguments can be scalar or matrix of the same size. Example: [E,N] = meshgrid(linspace(-10,10,50)); dG = okubo92(E,N,6,90,90,10,10,0,5,0,2670); figure, pcolor(E,N,1e8*dG), shading interp hold on, [c,h]=contour(E,N,1e8*dG,-50:10:50,'k'); clabel(c,h), hold off colorbar, polarmap reproduces the figure 4a in Okubo [1992] paper, considering a 10x10km fault at 6km depth, oriented N90°-strike, 90°-dip, and 5m senestral dislocation. Gravity changes are computed on a regular grid from -10 to 10km, and are plotted as a surface in 1e-8 m/s^2 unit. POLARMAP is an author's colormap. See also OKADA85 author's function to compute complete displacements, tilt and strain at free surface. Author: Francois Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2012-06-11 Updated: 2012-06-12 References: Aki K., and P. G. Richards, Quantitative seismology, Freemann & Co, New York, 1980. Okubo S., Gravity and Potential Changes due to Shear and Tensile Faults in a Half-Space, J. Geophys. Res., 97:B5, 7137-7144, 1992. Download okubo92.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
|
The Standard for the Exchange of Earthquake Data (SEED) is an international standard format for the exchange of digital seismological data. SEED was designed for use by the earthquake research community, primarily for the exchange between institutions of unprocessed earth motion data. It is a format for digital data measured at one point in space and at equal intervals of time. The SEED format consists of Volume Control Headers, Abbreviation Control Headers, Station Control Headers, Time Span Control Headers and finally Data Records. In complement to “Dataless” SEED volumes, exists the “Data-only” volume called Mini-SEED (see http://www.iris.edu for further information). The purpose of these functions is to read and write miniSEED data files directly from Matlab, avoiding intermediate format conversion (like SAC or other formats for which many functions exist), having a full control on headers and formats. Note that IRIS proposes a complete C-library for miniSEED (libmseed) and there is a free Matlab toolbox using this library to read and write miniSEED files (see Stefan Mertl site), but it needs some compilation to be used (unsignificant for geeks, quite easy on UNIX/Linux, but no doubt more problematic for Windows users...). The present functions rdmseed.m and mkmseed.m are 100% Matlab code standalone scripts, slower, but fully portable and platform independent. rdmseed: reading miniSEED file Each data record is imported into a structure array, allowing to adress data blocks and header fields individually (useful for multi-channel files), just as concatenating all data with a simple cat(1,X.d) function. Time stamps are also converted into Matlab datenum format. The function reads miniSEED "data-only" using the two mostly used compression formats Steim-1 and Steim-2. General FDSN formats have also been implemented (ASCII, 16/24/32-bit integers, IEEE floats and doubles), and GEOSCOPE multiplexed old formats (24-bit, 16/3 or 16/4-bit gain ranged). All these formats should work but some of them have not been tested using real data. I also partly coded Steim-3 format but without a clear description and any file example... Since I never met any data file using this format, I don't know if it's really useful. The function detects also automatically big/little-endian coded files. Known Blockettes are 1000, 1001, 100, 500 and 2000. If there is no Blockette 1000 (which is mandatory in SEED format...), default 4096-byte record length, big-endian and Steim-1 compression are used. These values can be set using additional arguments. Using extra output argument, some analysis can be done on the data stream (detection of gaps and overlaps), and channel components are detected. Without any output arguments, or with an additionnal 'plot' input argument, the function plots the imported data in a new figure (works also in case of multi-channel file). Steim-1/2 compression decoding strategy has been deeply optimized for Matlab. The proposed method, as vectorized as possible, is about 30 times faster than a 'C-like' loops coding... which is still 10 times slower than the same C-compiled program, but, well, this is the Matlab's other side of the coin! mkmseed: writing miniSEED file The function allows to export a data vector D to miniSEED file, giving origin date and time T0 and sampling rate FS (for monotonic data) or a time vector T. Header information is specified using the filename string with conventional naming "Network.Station.Location.Channel". Output file names will have appended ".Year.Day" and multiple file may be produced if data exceed a day. Data encoding format can be specified (16/32-bit integers, IEEE float/double, Steim-1/2; Geoscope 16/3-4). If not, it will depend on the class of variable D. Note that D will be converted to relevant class depending on encoding format. Be aware that Steim-2 compression uses 30-bit signed integer to code difference between 2 consecutive samples; equivalent to a 29-bit coding. Steim-1 uses 32-bit signed integers to code differences, so a 31-bit equivalent. This means that full range of signed 32-bit integer data cannot be correctly encoded with these formats in some exceptional cases (may lead to data loss). Finally, Steim-1/2 encoding needs time processing to achieve a maximum compression ratio of about 6 in file size, while other formats are very fast but result in large files. Binary file will be big-endian coded, with blockettes 1000, and default record length is 4096 bytes (this may be changed using input argument). Type "doc rdmseed" or "doc mkmseed" for detailed usage. MKMSEED Write data in miniSEED file format. MKMSEED(FILENAME,D,T0,FS) writes miniSEED file FILENAME from stritl monotonic data vector D, time origin T0 (a scalar in Matlab datenum compatible format) and sampling rate FS (in Hz). Encoding format will depend on D variable class (see below). MKMSEED(FILENAME,D,T,FS) where T is a time vector of the same length as data vector D, will create data records of monotonic blocks of samples, splitting each time the sampling frequency FS is not equal to time difference between two successive samples (with a 50% tolerency). Network, Station, Channel and Location codes will be extracted from FILENAME which must respect the format "NN.SSSSS.LC.CCC" where: NN = Network Code (2 characters max, see FDSN list) SSSSS = Station identifier (5 char max) LC = Location Code (2 char max) CCC = Channel identifier (3 char max) Final filename will have appended string ".YYYY.DDD" corresponding to year and ordinal day of origin time (from T0 value). Multiple files may be created if time span of data exceeds day limit. MKMSEED(...,EF,RL) specifies encoding format EF and data record length RL (in bytes). RL must be a power of 2 greater or equal to 256. Data encoding format EF must be one of the following FDSN codes: 1 = 16-bit integer (default for class 2-bit, 8-bit, int16) 3 = 32-bit integer (default for class uint16, int32) 4 = IEEE float32 (default for class single) 5 = IEEE float64 (default for all other class) 10 = Steim-1 compression (D will be converted to int32) 11 = Steim-2 compression (D will be converted to int32) 13 = Geoscope 16-3 gain ranged (D will be converted to double) 14 = Geoscope 16-4 gain ranged (D will be converted to double) MKMSEED(FILENAME,D) uses default value for T0 = now (present date and time), and FS = 100 Hz. File(s) will be coded big-endian, flags set to zero, blockette 1000, default record length is 4096 bytes. Outputs have been tested with PQL II software from IRIS PASSCAL (v2010.268). See also RDMSEED function for miniSEED file reading. Acknowledgments: Florent Brenguier, Julien Vergoz, Constanza Pardo. References: IRIS (2010), SEED Reference Manual: SEED Format Version 2.4, May 2010, IFDSN/IRIS/USGS, http://www.iris.edu IRIS (2010), PQL II Quick Look trace viewing application, PASSCAL Instrument Center, http://www.passcal.nmt.edu/ Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2011-10-19 Updated: 2014-05-07 : new version of rdmseed on April 21, 2012 with bug corrections (little-endian coding). Update your file! Download rdmseed.m and mkmseed.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). |
|
RDSAC reads a seismic data file encoded in the IRIS/SAC binary format, and returns a time vector, a data vector and all header variables in a structure (field names correspond to exact IRIS variable names). Type "doc rdsac" for detailed usage. RDSAC Read SAC data file. X=RDSAC(FILE) reads the Seismic Analysis Code (SAC) FILE and returns a structure X containing the following fields: t: time vector (DATENUM format) d: data vector (double) HEADER: header sub-structure (as defined in the IRIS/SAC format). RDSAC without input argument will open a file browser window. X=RDSAC(...,'plot') or RDSAC(...) without output argument will plot the data in a new figure. Reference: http://www.iris.edu/files/sac-manual/ Author: F. Beauducel <beauducel@ipgp.fr> Created: 2014-04-01 Updated: 2014-04-25 Download rdsac.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). |
|
|
This function plots regular grids of elevation X,Y,Z in a more efficient manner than SURFL Matlab's function, because it recomputes lighting and displays result as shaded color flat RGB image. It uses also median-style filter to deal with min/max values of elevations and gradient, and proposes two specific colormaps "landcolor" and "seacolor". Color mapping and lighting parameters can be changed from default values. In addition, several options are available: 'cartesian' to add decimal axis, 'latlon' to add geographical axis (GMT-like), 'legend' for an automatic scaling legend, 'lake' for automatic flat area color-filling and 'interp' to fill the novalue gaps... This may be useful to produce high-quality and moderate-size Postscript image adapted for publication. Figure examples:
See "doc dem" for syntax, examples and help. See also the READHGT function from same author. DEM Shaded relief image plot DEM(X,Y,Z) plots the Digital Elevation Model defined by X and Y coordinate vectors and elevation matrix Z, as a lighted image using specific "landcolor" and "seacolor" colormaps. DEM uses IMAGESC function which is much faster than SURFL when dealing with large high-resolution DEM. It produces also high-quality and moderate-size Postscript image adapted for publication. [H,I] = DEM(...); returns graphic handle H and illuminated image as I, an MxNx3 matrix (if Z is MxN and DECIM is 1). DEM(X,Y,Z,'Param1',Value1,'Param2',Value2,...) specifies options or parameter/value couple (case insensitive): --- Lighting options --- 'Azimuth',A Light azimuth in degrees clockwise relative to North. Default is A = -45 for a natural northwestern illumination. 'Contrast',C Light contrast, as the exponent of the gradient value: C = 1 for linear contrast (default), C = 0 to remove lighting, C = 0.5 for moderate lighting, C = 2 or more for strong contrast. 'LCut',LC Lighting scale saturation cut with a median-style filter in % of elements, such as LC% of maximum gradient values are ignored: LC = 0.2 is default, LC = 0 for full scale gradient. 'km' Stands that X and Y coordinates are in km instead of m (default). This allows correct lighting. Ignored if LATLON option is used. --- Elevation colorscale options --- 'ZLim',[ZMIN,ZMAX] Fixes min and max elevation values for colormap. Use NaN to keep real min and/or max data values. 'ZCut',ZC Median-style filter to cut extremes values of Z (in % of elements), such that ZC% of most min/max elevation values are ignored in the colormap application: ZC = 0.5 is default, ZC = 0 for full scale. --- "No Value" elevation options --- 'NoValue',NOVALUE Defines the values that will be replaced by NaN. Note that values equal to minimum of Z class are automatically detected as NaN (e.g., -32768 for int16 class). 'NaNColor',[R,G,B] Sets the RGB color for NaN/NoValue pixels (default is a dark gray). Note that your must specify a valid 3-scalar vector (between 0 and 1); color characters like 'w' or 'k' are not allowed, use [1,1,1] or [0,0,0] instead. 'Interp' Interpolates linearly all NaN values (fills the gaps using linear triangulation), using an optimized algorithm. --- Colormap options --- 'LandColor',LMAP Uses LMAP colormap instead of default (landcolor, if exists or jet) for Z > 0 elevations. 'SeaColor',SMAP Sets the colormap used for Z ≤ 0 elevations. Default is seacolor (if exists) or single color [0.7,0.9,1] (a light cyan) to simulate sea color. 'ColorMap',CMAP Uses CMAP colormap for full range of elevations, instead of default land/sea. This option overwrites LANDCOLOR/SEACOLOR options. 'Lake' Detects automaticaly flat areas different from sea level (non-zero elevations) and colors them as lake surfaces. 'Watermark',N Makes the whole image lighter by a factor of N. --- Basemap and scale options --- 'Cartesian' Plots classic basemap-style axis, considering coordinates X and Y as cartesian in meters. Use parameter "km' for X/Y in km. 'LatLon' Plots geographic basemap-style axis in deg/min/sec, considering coordinates X as longitude and Y as latitude. Axis aspect ratio will be adjusted to approximatively preserve distances (this is not a real projection!). This overwrites ZRatio option. 'Legend' Adds legends to the right of graph: elevation scale (colorbar) and a distance scale (in km). 'FontSize',FS Font size for X and Y tick labels. Default is FS = 10. 'BorderWidth',BW Border width of the basemap axis, in % of axis height. Must be used together with CARTESIAN or LATLON options. Default is BW = 1%. --- Decimation options --- For optimization purpose, DEM will automatically decimate data to limit to a total of 1500x1500 pixels images. To avoid this, use following options, but be aware that large grids may require huge computer ressources or induce disk swap or memory errors. 'Decim',N Decimates matrix Z at 1/N times of the original sampling. 'NoDecim' Forces full resolution of Z, no decimation. --- Informations --- Colormaps are Mx3 RGB matrix so it is easy to modify saturation (CMAP.^N), set darker (CMAP/N), lighter (1 - 1/N + CMAP/N), inverse it (flipud(CMAP)), etc... To get free worldwide topographic data (SRTM), see READHGT function. For backward compatibility, the former syntax is still accepted: DEM(X,Y,Z,OPT,CMAP,NOVALUE,SEACOLOR) where OPT = [A,C,LC,ZMIN,ZMAX,ZC], also option aliases DEC, DMS and SCALE, but there is no argument checking. Please prefer the param/value syntax. Author: François Beauducel <beauducel@ipgp.fr> Created: 2007-05-17 Updated: 2013-03-10 Download dem.m, landcolor.m and seacolor (needed) at Matlab Central file exchange (Copyright 2013 F. Beauducel / BSD License). |
|
|
This function imports .HGT "height" binary data files from NASA SRTM global digital elevation model of Earth land, corresponding to 1x1 degree tiles of 3-arc seconds resolution (SRTM3, around 90 m) and 1-arc second (SRTM1, around 30 m) for USA territory, and returns coordinates vectors latitude and longitude, and a matrix of elevation values. The function includes an automatic download of data from the USGS SRTM webserver, so indicating latitude and longitude is sufficient to get the data and instant map plot anywhere in the World. Also some basic options as 'merge' to concatenate tiles, 'interp' for gaps (novalue) linear interpolation, 'crop' for rectangle selection inside tile(s). READHGT Import/download NASA SRTM data files (.HGT). X=READHGT(FILENAME) reads HGT "height" SRTM data file FILENAME and returns a structure X containing following fields lat: coordinate vector of latitudes (decimal degree) lon: coordinate vector of longitudes (decimal degree) z: matrix of elevations (meters, INT16 class) FILENAME must be in the form "[N|S]yy[E|W]xxx.hgt[.zip]" as downloaded from SRTM data servers. X=READHGT(LAT,LON) attemps to download the file *.hgt.zip corresponding to LAT and LON (in decimal degrees) coordinates (lower-left corner) from the USGS data server (needs an Internet connection and a companion file "readhgt_srtm_index.txt"). Downloaded filename(s) will be given in an additional output structure field X.hgt. LAT and/or LON can be vectors: in that case, tiles corresponding to all possible combinations of LAT and LON values will be downloaded, and optional output structure X will have as much elements as tiles. READGHT(...) without output argument or X=READHGT(...,'plot') plots the tile(s). For better plot results, it is recommended to install DEM personal function available at author's Matlab page. READHGT(LAT,LON,...,'merge'), in case of adjoining values of LAT and LON, will concatenate tiles to produce a single one. ATTENTION: tiles assembling may require huge computer ressources and cause disk swaping or memory error. READHGT(LAT,LON,...,'interp') linearly interpolates missing data. READHGT(...,'crop',[LAT1,lAT2,LON1,LON2]) crops the map using latitude/longitude limits. READHGT(LAT,LON,...,'crop'), without limits argument vector, crops the resulting map around existing land (reduces any sea or novalue areas at the borders). READHGT(LAT,LON,...,'srtm3') forces SRTM3 download (by default, SRTM1 tile is downloaded if exists). Usefull for USA neighborhood. READHGT(LAT,LON,OUTDIR) specifies output directory OUTDIR to write downloaded files. READHGT(LAT,LON,OUTDIR,URL) specifies the URL address to find HGT files (default is USGS). Examples: - to plot a map of the Paris region, France (single tile): readhgt(48,2) - to plot a map of Flores volcanic island, Indonesia (5 tiles): readhgt(-9,119:123,'merge') - to download SRTM1 data of Cascade Range (27 individual tiles): X=readhgt(40:48,-123:-121); Informations: - each file corresponds to a tile of 1x1 degree of a square grid 1201x1201 of elevation values (SRTM3 = 3 arc-seconds), and for USA territory at higher resolution 3601x3601 grid (SRTM1 = 1 arc-second) - elevations are of class INT16: sea level values are 0, unknown values equal -32768 (there is no NaN for INT class), use 'interp' option to fill the gaps - note that borders are included in each tile, so to concatenate tiles you must remove one row/column in the corresponding direction (this is made automatically with the 'merge' option) - downloaded file is written in the current directory or optional OUTDIR directory, and it remains there - NASA Shuttle Radar Topography Mission [February 11 to 22, 2000] produced a near-global covering on Earth land, but still limited to latitudes from 60S to 60N. Offshore tiles will be output as flat 0 value grid - if you look for other global topographic data, take a look to ASTER GDEM, worldwide 1 arc-second resolution (from 83S to 83N): http://gdem.ersdac.jspacesystems.or.jp (free registration required) Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris References: http://dds.cr.usgs.gov/srtm/version2_1 Acknowledgments: Yves Gaudemer Created: 2012-04-22 Updated: 2013-01-17 Download readhgt.m and readhgt_srtm_index.txt (needed) at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
|
UTM2LL converts Universal Transverse Mercator (UTM) East/North coordinates to latitude/longitude. LL2UTM converts latitude/longitude coordinates to UTM. Both functions are using precise formula (millimeter precision), possible user-defined datum (WGS84 is the default), and are all vectorized (no loop in the code). It means that huge matrix of points, like an entire DEM grid, can be converted very fast. Example (needs readhgt.m author's function):
X = readhgt(36:38,12:15,'merge','crop',[36.5,38.5,12.2,16],'plot');
[lon,lat] = meshgrid(X.lon,X.lat);
[x,y,zone] = ll2utm(lat,lon);
z = double(X.z);
z(z==-32768 | z<0) = NaN;
figure
pcolor(x,y,z); shading flat; hold on
contour(x,y,z,[0,0],'w')
hold off; axis equal; axis tight
xlabel('East (m)'); ylabel('North (m)')
title(sprintf('Sicily - UTM zone %d WGS84',zone))
loads SRTM full resolution DEM of Sicily in lat/lon (a 2400x4500 grid), converts it to UTM and plots the result. To make a regular UTM grid, you may interpolate x and y with griddata function. See "doc ll2utm" and "doc utm2ll" for syntax and help. LL2UTM Lat/Lon to UTM coordinates precise conversion. [X,Y]=LL2UTM2(LAT,LON) converts coordinates LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84. LAT and LON can be scalars, vectors or matrix. Outputs X and Y will have the same size as inputs. LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be a string in the following list: 'wgs84': World Geodetic System 1984 (default) 'nad27': North American Datum 1927 'clk66': Clarke 1866 'nad83': North American Datum 1983 'grs80': Geodetic Reference System 1980 'int24': International 1924 / Hayford 1909 or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in meters) and F is flattening of the user-defined ellipsoid. [X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value stands for southern hemisphere points). Notice: - LL2UTM does not perform cross-datum conversion. - precision is near a millimeter. Reference: I.G.N., Projection cartographique Mercator Transverse: Algorithmes, Notes Techniques NT/G 76, janvier 1995. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2003-12-02 Updated: 2014-02-17 UTM2LL UTM to Lat/Lon coordinates precise conversion. [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters) defined in the UTM ZONE (integer) to latitude LAT and longitude LON (in degrees). Default datum is WGS84. X and Y can be scalars, vectors or matrix. Outputs LAT and LON will have the same size as inputs. For southern hemisphere points, use negative zone -ZONE. UTM2LL(X,Y,ZONE,DATUM) uses specific DATUM for conversion. DATUM can be a string in the following list: 'wgs84': World Geodetic System 1984 (default) 'nad27': North American Datum 1927 'clk66': Clarke 1866 'nad83': North American Datum 1983 'grs80': Geodetic Reference System 1980 'int24': International 1924 / Hayford 1909 or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in meters) and F is flattening of the user-defined ellipsoid. Notice: - UTM2LL does not perform cross-datum conversion. - precision is near a millimeter. Reference: I.G.N., Projection cartographique Mercator Transverse: Algorithmes, Notes Techniques NT/G 76, janvier 1995. Author: Francois Beauducel, <beauducel@ipgp.fr> Created: 2001-08-23 Updated: 2014-02-17 Download ll2utm.m and utm2ll.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). |
 
|
This program computes a map of visibility from a selected point on a topography. It has been written to help the search for radio relay best location. Because it considers only direct line of sight, it gives a good estimation for possible radio link for short distances only (less than 10 km), neglecting curvature of the Earth, Fresnel zone and atmospheric refraction on radio waves propagation. The program computes the relative elevation angle of the mask for each point (the angle is null or negative if the point is visible). The function needs a digital elevation model Z and associated (X,Y) vectors or matrices of coordinates (same unit as Z), position of the point (X0,Y0), the antenna height H0 (for instance 4 m), and the hypothetic antenna height Ha on each topography points (for instance 3 m). When no output argument is given, the function plots a map of the results (color map of mask angles, and blank for visible points, see example screenshot). The script is not fully optimized because it makes a global loop on the matrix elements to compute each profile of topography... (I didn't find (yet) a way to fully vectorize the problem), so it takes some time to compute, depending on the number of element of Z. However, I found a faster algorithm (about 2 times faster), giving approximate result but usefull to process a first map. See help for syntax, and script comments for details. RADIOCOVER Radio link coverage map. RADIOCOVER(X,Y,Z,X0,Y0,H0,Ha) computes the coverage map of possible direct linear radio link from the point (X0,Y0) with antenna height H0, using digital terrain model defined by coordinate vectors X and Y, and elevation matrix Z, and hypothetic antenna height Ha, then plots a color map of the relative elevation mask angle (in degrees) with blank areas where there is no mask (visible), together with a contour map of the topography. X and Y can be vectors with length(X) = n and length(Y) = m where [m,n] = size(Z), or matrices of the same size as Z (as from MESHGRID). RADIOCOVER(...,METHOD) specifies alternate methods. The default is nearest neighbor interpolation. Available methods are: 'nearest' - nearest neighbor interpolation (default) 'linear' - linear interpolation (smoother result) 'fast' - approximate algorithm (about 2 times faster) M = RADIOCOVER(...); returns a matrix of relative elevation mask angle (in degrees, same size as Z), without producing graphic. Visible points have null or negative values. The model assumes linear propagation of radio waves (direct line of sight between the two antennas), and neglects curvature of the Earth, Fresnel zone, and atmospheric refraction. Example: [x,y,z]=peaks(100); [fx,fy]=gradient(z); z=sqrt(fx.^2+fy.^2); surf(x,y,z), shading flat, light, view(-24,74) radiocover(x,y,z,-0.84,-0.27,.05,.05) Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2003-01-10 Modified: 2010-01-08 Download radiocover.m at Matlab Central file exchange (Copyright 2003-2009 F. Beauducel / BSD License). |
 |
The function GREATCIRCLE computes the shortest path along the great circle ("as the crow flies") between two points defined by their geographic coordinates (latitude and longitude). With one output argument it returns distance or vector of distances, with two or more output arguments it returns path coordinates and optional vector of distances and bearing angles along the path. The function LOXODROME computes the path with a constant bearing, crossing all meridians of longitude at the same angle. It returns also a vector of distances and the bearing angle. Loxodrome path (also known as "rhumb line") is longer than great circle one, but still used in navigation as it is easier to follow with a compass. Type 'doc greatcircle' or 'doc loxodrome' for syntax, help and examples. Download greatcircle.m and loxodrome.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
 |
This little function draws a compass rose with optional Cardinal directions with basic properties:
Download comprose.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
 |
This little function could be an alternative to other existing arrow plot functions, since it has very simple, vectorized and effective coding. In fact, arrows are all plotted in a single patch command! The function also accepts any standard patch parameter/value pair like 'FaceColor', 'EdgeColor', 'LineWidth', etc. It can be used to plot a personalized arrow at coordinates X,Y with length L and azimuth AZ, or any of these arguments can be vector or matrix, like QUIVER function, to plot multiple arrows. Arrow's aspect ratio, head and segment line shapes are configurable with 4 parameters: head width, head length, head inside length and segment line width, all normalized to arrow's length. Type 'doc arrows' for syntax, help and examples. Download arrows.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
|
|
This function gives the parameters for about 40 main wave components of Earth tides due to Sun and Moon (annual, semiannual, monthly, fortnightly, diurnal, semidiurnal, ...) using the harmonic development of Arthur Thomas Doodson [1890-1968]. For a given wave, it returns in a structure the 6 Doodson's arguments, wave Darwin's symbol, wave description and the wave period (in days).
>> X=doodson('M2')
X =
symbol: {'M2'}
name: {'Principal lunar semidiurnal'}
doodson: [2 0 0 0 0 0]
period: 0.5175
Type doodson without argument to see the list of available waves.
Download doodson.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). |
|
|
This function reads data files from Fluke TI32 infrared thermal camera (.IS2 proprietary format). It returns data images (raw and calibrated temperatures, and visible RGB image) and some of the header information in a single structure. Without output arguments, it plots blended infrared and visible images and a histogram of temperatures. The function has been tested for TI32 camera original files only. Presently it does not work with files edited with the Fluke software, neither for other cameras. Also the function needs a separated calibration file since I was not able to find the original tables in the binary file (but I'm quite sure they exist...). Any help to improve that is welcome. READIS2 Import data files from Fluke infrared thermal camera (IS2 format). X=READIS2(FILENAME) reads IS2 radiometric file FILENAME and returns a structure X containing fields with data and headers: IR: infrared sensor values (240x320 uint16 matrix) TB: brillancy temperatures, in °C (240x320 double matrix) T: body temperatures, in °C (240x320 double matrix) VI: visible image (480x640x3 uint8 RGB) VI2: cropped and interpolated visible image adjusted on IR image pixels (240x320x3 uint8 RGB) Emissivity: emissivity value (0 to 1.00) Transmission: transmission value (0 to 1.00) BackgroundTemp: background temperature (°C) Date: [YEAR MONTH DAY HOUR MINUTE SECOND] X=READIS2(FILENAME,MAP) plots temperature image using MAP color scale over visible image. MAP is a valid [n x 3] RGB colormap, for instance JET or GRAY. Use MAP=[] to specify other parameters without plot. X=READIS2(FILENAME,MAP,PARAM) applies given parameters to overwrite values from camera. PARAM = [DT,EM,TR,BT] (use NaN to keep original values): DT = time difference (in days) applied to image timestamp EM = emissivity TR = transmission BT = background temperature (°C) X=READIS2(FILENAME,MAP,PARAM,CAL) applies a polynomial function to calibrate brillancy temperatures, as follow: TB_cal = polyval(CAL,TB_raw). This will affect also body temperatures. Example: CAL = [1.08,-2.72] is for my own Ti32 camera calibration result for 0-100°C range laboratory experiment. Notes: - the function has been tested only with original IS2 files from a Fluke TI32 camera. - to produce temperatures from IR sensor values, the function needs a calibration file 'readis2_ti32.dat' located in current directory or Matlab path (see corresponding header info). - body temperatures Tr (output X.T) and brillancy temperatures Tb (output X.IR) are related with the formula: Tb = ε*τ*Tr + (2-ε-τ)*Te where ε is emissivity, τ is transmission, and Te is background (reflected) temperature. So Tb and Tr are equal for ε = τ = 1. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2011-07-24 Updated: 2012-04-01 Acknowledgments: Magnus Andersen (read_is2.m), Pierre Agrinier, Valérie Clouard, Damien Gaudin, David Krammer. Download readis2.m and readis2_ti32.dat (needed) here (Copyright 2012 F. Beauducel / BSD License). |
|
When representing "polarized" data for which zero and sign are meaningful (like for instance, a gravity anomalies map), it is useful to have a colormap where white stands for zero value, and greater magnitudes are plotted as colors with increasing saturation. This little function POLARMAP has two different usages:
POLARMAP Polarized color map POLARMAP applies a "polarized" red-white-blue colormap to current figure, and adjusts the color axis limits to be centered to zero. POLARMAP(M) fixes the number of colors to M (default is 64). POLARMAP(MAP) applies linear shading to white to the center of colormap MAP which can be any of existing colormaps (an Mx3 matrix of RGB). POLARMAP(MAP,C) uses exponent C to modify the shading contrast. Default is C = 1 for linear shading. Use C = 2 to strengthen the shading, or C = 0.5 to attenuate it. C=POLARMAP(...) returns an M-by-3 matrix containing the colormap, that can be used with COLORMAP function like other colormaps. Examples: pcolor(peaks), shading interp polarmap, colorbar then try the following polarmap(jet,0.5) Note the polar shading has no real interest with colormaps that include white color as one of the extremes (like GRAY, BONE, HOT, ...). See also JET, HSV, COPPER, SPRING, SUMMER, WINTER, COOL, COLORMAP, RGBPLOT. Author: François Beauducel <beauducel@ipgp.fr> Created: 2011-10-26 Updated: 2012-06-09 Download polarmap.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
|
This function is a simple digital filter with one's coefficients, nearly equivalent to moving average but instead of mean values, it computes the sum. This is particularly usefull for rainfall plots for which it is, in my opinion, the only good way to show the data. The figure attempts to demonstrate it. Let us have some rainfall data at 10-minute sampling interval, and simply want to plot the 1-hour rainfall... If you use classical histogram representation, it leads to artifacts that depends on a priori time windows offset you chose (see figure):
MOVSUM Moving sum plot. Y = MOVSUM(X,N) returns the "moving sum" of X on N consecutive values. Y has the same size as X, and each element of Y is the sum of the N previous values of X. If X is a regularly spaced data vector, then Y corresponds to a digital filter with N b-coefficient equal to one: y(n) = x(n) + x(n-1) + ... + x(1) This is nearly equivalent to a moving average, but instead of mean it is the sum. If X is a matrix, MOVSUM works down the columns. A typical use is for rainfall data vectors. For instance, if X contains rainfall values at 1-minute interval, MOVSUM(X,60) will be the hourly rainfall, expressed as a continous vector of the same size as X. Compared to classical histogram representation, this method avoids some artifacts due to windowing using a priori time intervals. MOVSUM(...) without output arguments produces a plot of the results. Author: François Beauducel <beauducel@ipgp.fr> Created: 2005 Updated: 2010-10-10 Download movsum.m (Copyright 2010 F. Beauducel / BSD License). |
|
This little function rounds a number (or the elements of a vector or matrix) towards the nearest number with N significant digits. This is a useful complement to Matlab's ROUND, ROUND10 and ROUNDN (Mapping toolbox), especially when dealing with data with a large variety of order of magnitudes. ROUNDSD Round with fixed significant digits ROUNDSD(X,N) rounds the elements of X towards the nearest number with N significant digits. ROUNDSD(X,N,METHOD) uses following methods for rounding: 'round' - nearest (default) 'floor' - towards minus infinity 'ceil' - towards infinity 'fix' - towards zero Examples: roundsd(0.012345,3) returns 0.0123 roundsd(12345,2) returns 12000 roundsd(12.345,4,'ceil') returns 12.35 See also Matlab's functions ROUND, ROUND10, FLOOR, CEIL, FIX, and ROUNDN (Mapping Toolbox). Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Acknowledgments: Edward Zechmann, Daniel Armyr, Yuri Kotliarov Created: 2009-01-16 Updated: 2015-04-03 Download roundsd.m at Matlab Central file exchange (Copyright 2015 F. Beauducel / BSD License). NEWS: This file was selected as MATLAB Central Pick of the Week, and integrated as an option of the core function ROUND in the new Matlab 2014b ! |
|
A little function that I imagined many years ago: a generalized "median-style" filter, an elegant way to generalize MIN, MAX, MEDIAN and extreme filtering functions, reduced to a single parameter. Let's sorting all values of a vector X in ascending order, thus consider the normalized rank position N of elements:
A typical example of use is minmax(X,[0.01 0.99]) which returns minimum and maximum values of X(:) but excluding the 1% amount of extreme values. This is particularily efficient for automatic scaling of noisy data (see the screenshot example), compared to the use of MEAN and STD functions which is biased by any high-magnitude values in X. Default behavior of MINMAX is to return a vector of minimum and maximum values. So minmax(X) is an abreviation of minmax(X,[0 1]), and equivalent of [min(X(:)) max(X(:))]. Compared to the built-in functions MIN and MAX, it does not return a vector from matrix but returns a scalar by considering all elements of matrix X. Also, it removes any NaN values (returns NaN only if all X values are NaN). Type 'doc minmax' for syntax, help and other examples. MINMAX Generalized median filter. Y=MINMAX(X) returns a 2 element vector Y = [min(X) max(X)]. Minimum and maximum values of X(:) are computed after excluding any NaN values. Y=MINMAX(X,N) with N between 0 and 1, returns generalized median filter value of X(:), e.g., it sorts elements of X then interpolates at (100*N) % rank position. N can be scalar, vector or matrix, so Y will have the same size as N. Examples: MINMAX(X,0) returns the minimum value of X elements. MINMAX(X,1) returns the maximum value of X elements. MINMAX(X,0.5) returns the median value of X elements. This is the equivalent of MEDIAN(X(:)). MINMAX(X,0.9) returns maximum value of X after excluding the 10% highest value elements. MINMAX(X,[0 1]) is the same as MINMAX(X). MINMAX(X,[0.01 0.99]) is a convenient way to compute automatic scale of X when samples are noisy, since it filters the 1% elements with extreme values. It may be used for color scaling with CAXIS or for plot scaling with set(gca,'YLim',...). See also MIN, MAX and MEDIAN. Author: François Beauducel <beauducel@ipgp.fr> Created: 1996 Updated: 2012-10-10 Download minmax.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). |
|
|
This is the core formula of discrete Fourier transform: it simply computes the amplitude and phase shift of fundamental or harmonics of a phase signal. Example: t = linspace(0,2*pi); x = 2*cos(t + pi/2) - cos(3*t) + rand(size(t)); harmfit(t,x,1:4)returns estimations of amplitudes/phases for the first four harmonics. HARMFIT Sinusoidal harmonic curve fitting. H = HARMFIT(X,Y) returns first harmonic amplitude and phase of the data vector Y relative to phase vector X (in radian) in a 3-element vector H = [N,AMP,PHA], so that AMP*COS(N*X + PHA) fits Y. HARMFIT(X,Y,N) computes the N'th harmonic. N = 1 stands for fundamental, N = 2 is second harmonic, etc... N can be a scalar or a vector of positive integers. For example use N = 1:4 to compute the first four harmonics. [H,YY] = HARMFIT(...) returns also an evaluation of harmonic curve fit in vector YY (as the sum of harmonics in N). This is simply the core calculation of discrete Fourier transform. Example: t = linspace(0,2*pi); x = 2*cos(t + pi/2) - cos(3*t) + rand(size(t)); harmfit(t,x,1:4) returns estimations of amplitudes/phases for the first four harmonics. Note that negative amplitudes are fitted with positive value with a PI phase difference. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2014-05-22 Updated: 2014-05-24 Download harmfit.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). |
|
These two scripts convert Roman numerals to and from any integer (scalar, vector or matrix), including large numbers greater than 4999 with the parenthesis notation (multiplies by 1000). The function NUM2ROMAN uses strict rules of modern notation (substractive principle for 4 and 9 bases) except for the common 'MMMM' form replacing '(IV)'. ROMAN2NUM is more flexible and is able to convert some other Roman notation possibilities, for instance the 3 different expressions roman2num({'IC','XCIX','XCVIIII'}) return [99,99,99], or roman2num({'MDXV','MCCCCCXV'}) return [1515,1515]. See help for syntax, and script comments for details. NUM2ROMAN Roman numerals. NUM2ROMAN(N) returns modern Roman numeral form of integer N, which can be scalar (returns a string), vector or matrix (returns a cell array of strings, same size as N). The function uses strict rules whith substractive notation and commonly found 'MMMM' form for 4000. It includes also parenthesis notation for large numbers (multiplication by 1000). Examples: num2roman(1968) num2roman(10.^(0:7)) reshape(num2roman(1:100),10,10) See also ROMAN2NUM. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2005 Modified: 2009-12-21 Download num2roman.m and roman2num.m at Matlab Central file exchange (Copyright 2005-2009 F. Beauducel / BSD License). |
|
This little function computes the date of Easter Sunday, using Oudin's algorithm for Julian calendar (starting 325 AD) and Gregorian calendar (after 1583 AD). Easter day is usefull to calculate other christian feasts, like Ash Wednesday (Easter - 46), Good Friday (Easter - 2), Ascension Thursday (Easter + 39), Pentecost (Easter + 49). Some translations of Easter: Pâques (French), Pascua (Spanish), Ostern (German), Pasqua (Italian), Paskah (Indonesian), Wielkanoc (Polish) ... See help for syntax, and script comments for details. EASTER Easter Day EASTER displays the date of Easter Sunday for present year. EASTER(YEAR) displays the date of Easter for specific YEAR, which can be scalar or vector. E = EASTER(...) returns Easter day(s) in DATENUM format. This function computes Easter Day using the Oudin's algorithm [1940], which is valid for Catholic Easter day from 325 AD (beginning of the Julian calendar). Easter day is usefull to calculate other usual christian feasts: datestr(easter-47) is Mardi Gras datestr(easter-46) is Ash Wednesday datestr(easter-24) is Mi-Carême datestr(easter-2) is Good Friday datestr(easter+1) is Easter Monday datestr(easter+39) is Ascension Thursday datestr(easter+49) is Pentecost Reference: Oudin, 1940. Explanatory Supplement to the Astronomical Almanac, P. Kenneth Seidelmann, editor. Tondering, C, 2008. http://www.tondering.dk/claus/calendar.html Author: François Beauducel, <beauducel@ipgp.fr> Created: 2002-12-26 Updated: 2009-03-01 Download easter.m at Matlab Central file exchange (Copyright 2002-2009 F. Beauducel / BSD License). |
|
The current Javanese calendar was inaugurated by Sultan Agung of Mataram in the Gregorian year 1633. It is based on a combination of the Hindu calendar "Saka" and the Islamic calendar based on the lunar month, and contains different cycles: Pasaran (5-day), Dina Pitu (7-day), Wetonan (35-day), Mangsa (solar month), Wulan (Moon month), Pawukon (210-day), Tahun (Moon year), Windu (8-year), Kurup (120-year). Coincidences of these multiple cycles have special mystical meanings for any Javanese people, for instance the birthday "Weton" or the Noble Days "Dino Mulyo". This is the primary time-keeping system for all matters having cultural, historical, and metaphysical significance in the Java island, Indonesia. This little script computes dates in the Javanese calendar, indicating Dina Pitu, Pasaran, Dino, Wulan, Tahun, Windu, and Kurup for today or a specific list of days. It indicates also the Noble Days if necessary. If you specify your date of birth, it can calculate your "weton" and a list of your javanese birthdays. For example, weton(1968,12,3) returns: A second use of this function is to display a month calendar that combines the 5-day "Pasaran" cycle and the 7-day Gregorian/Islamic week, called "Wetonan". See help for syntax, and script comments for details. Download weton.m at Matlab Central file exchange (Copyright 1999-2011 F. Beauducel / BSD License). See also my little program weton.zip (Perl) that does almost the same job (and more). |
