Nothing
`get.corner` <-
function( INfreq, INspec, dt, f1, f2, PLOT=FALSE, VERBOSE = FALSE)
{
#### calculate the best fit Omega0, Corner freq and Tstar
## for a Brune Model
## computations are done in the LOG-LOG domain
## so we fit a flat line for Omega and a linear fall off for tstar
## example: mc = get.corner(Spec$f, lspec, dt, 0.01, 15.0, PLOT=TRUE)
if(missing(PLOT)) { PLOT=FALSE }
if(missing(VERBOSE)) { VERBOSE = FALSE}
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
n = length(INfreq);
corn = 0
ave = 0
slope = 0
interc = 0
K = 0
flag = INfreq>=f1 & INfreq<=f2
freqlim = INfreq[flag]
speclim = INspec[flag]
Lfreqlim = log10(freqlim)
Lspeclim = log10(speclim)
n = length(freqlim);
barf<-.C("CALL_DCORN", PACKAGE = "RSEIS",
as.double(Lfreqlim),
as.double(Lspeclim),
as.integer(n) ,
as.integer(K) ,
as.double(ave) ,
as.double(slope) ,
as.double(interc)
)
K = unlist(barf[[4]])+1
ave=unlist(barf[[5]]);
slope = unlist(barf[[6]]);
interc = unlist(barf[[7]]);
### now get the corner frequency
newx = (ave-interc)/slope
## fc = freqlim[K]
fc = 10^newx
PI = 3.14159265358979;
gamma = 2.0;
f = fc;
ftem = f*(PI * log10(2.718281828));
fa = slope*log10(fc);
fb = 0.5*log10( 1+ (f/fc)^(2*gamma));
tem = -(fa+fb)/ftem;
tstar0 = tem/10.;
omega0 = 10^ave;
if(VERBOSE)
{
message(paste(sep=' ', "Gcorn:",
formatC(ave) ,
formatC(omega0),
formatC(fc),
formatC(slope),
formatC(interc),
formatC(tstar0) ))
}
if(PLOT==TRUE)
{
opar = par("usr")
par(mfrow=c(2,1))
plot(Lfreqlim, Lspeclim, type='l', main="Restricted Freq", xlab="Log Freq", ylab="Log Displacement")
abline(h=ave, col=4)
abline(interc, slope, col=2)
plot(log10(INfreq), log10(INspec), type='l', main="Full Freq", xlab="Log Freq", ylab="Log Displacement")
abline(h=ave, col=4)
abline(interc, slope, col=2)
invisible( par(opar))
}
ret = list(corn=fc,ave=ave,slope=slope,interc=interc,tstar0=tstar0,omega0=omega0 )
return(ret)
## example: mc = get.corner(Spec$f, lspec, dt, 0.01, 15.0, PLOT=TRUE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.