Monte Carlo simulation of exoplanet transit light curves in R

I've posted about the Kepler mission before, but this is a different approach.  While it's easy to find all sorts of animations and diagrams showing how the light from a star is blocked over time as a planet transits in front of it, I got to thinking how those "ideal' light curves are found.  I figured it was probably possible with some sort of convolution, but I wasn't really sure.  Finding a function to determine the area of an intersection of two moving circles didn't sound like a lot of fun anyways.  But, I thought it would be interesting to attack it as a simulation.

Since I've been fooling around with R lately, that's what I used.  My code is below.  It is a very simplified model, simulating one flat disk transiting another.  The code produces random points in a box around the "planet" and discards the ones not within the planet's radius.  Those points are then moved to the appropriate place in the transit across the "star", as determined by the step in the loop.  The points are checked to see if they lay within the star's radius, and that count is divided by the total number of points to determine the proportion of the planet's area that is intersecting the star.  That area is then divided by the star's area to find how much of the star's light is blocked.  Those proportions are collected at each point of the planet's trip and can then a simulated light curve can be charted.

The ratio of planet to star radius can be changed, as can the crossing latitude.  The simulation parameters of number of points generated and transit step length can also be changed.  Again, this is nothing fancy and there are probably much better ways of doing this, but it was a fun project to chew on.

star.radius <- 1  # Conceptually, it's easier if this is 1
planet.radius <- .20

star.crossing.lat <- 88
lat.radians <- star.crossing.lat * ((2 * pi) / 360)
star.crossing.y <- star.radius * sin(lat.radians)

step.length <- .001
planet.x.pos <- seq((-2 * star.radius), (2 * star.radius), step.length)  # Position of center of planet position, in terms of stellar radius.
steps <- as.integer(length(planet.x.pos))

point.count <- 25000
box.point.count <- ceiling(point.count * (4 / pi) * 1.1)  # Can be wasteful

light.curve <- vector(length = steps, mode = "numeric")

counter <- 1
for(i in planet.x.pos){
       
  #Generate poins on planet
  box.points_df <- data.frame(runif(box.point.count, min = -planet.radius, max = planet.radius), runif(box.point.count, min = -planet.radius, max = planet.radius))
  colnames(box.points_df) <- c('x', 'y')
  box.points_df$r <- sqrt(box.points_df$x ^ 2 + box.points_df$y ^ 2)
  box.points_df$in.circle <- box.points_df$r <= planet.radius
  circle.points_df <- box.points_df[box.points_df$in.circle == TRUE, ]
  circle.points_df <- circle.points_df[1:point.count, ]
   
  #Move planet points into position
  circle.points_df$x <- circle.points_df$x + i
  circle.points_df$y <- circle.points_df$y + star.crossing.y

  #Check if points intersect with star
  planet.point.check <- sqrt((circle.points_df$x) ^ 2 + (circle.points_df$y) ^ 2)
  intersect.points <- sum(planet.point.check < star.radius)
 
  #Calculate relative flux
  pct.planet.intersect <- intersect.points / nrow(circle.points_df)
  area.intersect <- pct.planet.intersect * pi * planet.radius ^ 2
  relative.flux <- ((pi * star.radius ^ 2) - area.intersect) / (pi * star.radius ^ 2)
  light.curve[counter] <- relative.flux

  counter <- counter + 1
  }


And here is a sample of the output for a transit for the parameters in the code above.