##### This is for fitting 
    ##### lambda(t,x,y) = mu rho(x,y) + K SUM gt(t-t_i)gxy(x-x_i,y-y_i)gmi(m_i), 
    ##### with rho(x,y) = 1/(X1Y1), 
    ##### gt(t) = (p-1)c^(p-1) (t+c)^(-p),
    ##### g(x,y) = (q-1)d^(q-1)/pi (r^2+d)^(-q), with x^2+y^2=r^2,
    ##### and gmi(m) = exp(a(m-M0)).  
    ##### The space S = [0,X1] x [0,Y1], in time [0,T], and lower
    ##### magnitude cutoff M0.
    ##### The magnitude density is b exp(-b (m-M0)).    
    ##### For any theta, the integral of lambda over the space time region is approx. 
    ##### mu T + K SUM exp{a(m_i-M0)}.
    ##### We can estimate b, for the magnitude frequency, separately.
    ##### Here the parameter vector theta = (mu, K, c, p, a, d, q)

## If you have: 
## datax 
## datay
## datat 
## datan 
## datam, 
## then write  
## z = list()
## z$t = datat
## z$lon = datax 
## z$lat = datay
## z$m = datam  
## z$n = datan 
## Also define T,X1,Y1, and M0 externally before calling this. 

## r2 will be the matrix of squared distances between pts. 
## dt will be the matrix of time differences between pts. 
r2 = matrix(0,nrow=z$n, ncol=z$n); dt = r2
for(i in 1:z$n){ for(j in 1:z$n){ 
r2[i,j] = (z$lon[j]-z$lon[i])^2+(z$lat[j]-z$lat[i])^2
dt[i,j] = t[j]-t[i]}}

m3 = function(x) signif(x,3) 

## This is the loglikelihood function in R. 
logletas = function(theta){
mu = theta[1]; K = theta[2]; c = theta[3]; p = theta[4] 
a = theta[5]; d = theta[6]; q = theta[7] 
cat("\n mu = ",m3(mu),", K = ",m3(K),", c = ",m3(c),", p = ",m3(p), 
", a = ",m3(a),", d = ",m3(d),", q = ",m3(q),"\n") 
if(min(mu,K,c,p,a,d,q)<0.000000001) return(99999) 
if(K>.99999) return(99999)
sumlog = log(mu/X1/Y1) 
intlam = mu*T + K*sum(exp(a*(z$m-M0))) 
const = K*(p-1)*c^(p-1)*(q-1)*d^(q-1)/pi 
for(j in 2:(z$n)){
   lamj = mu/X1/Y1 + const*sum(
	exp(a*(z$m[1:(j-1)]-M0))*(dt[1:(j-1),j]+c)^(-p)*(r2[1:(j-1),j]+d)^(-q)) 
   if(lamj < 0){
    cat("lambda ",j," is less than 0.")
    return(99999)
   }
   sumlog = sumlog + log(lamj)
}
loglik = sumlog - intlam
cat("loglike is ", loglik, ". sumlog = ", sumlog,". integral = ", intlam,".\n")
return(-1.0*loglik)
}