library(FITSio)
library(data.table)
library(ggplot2)

Solar Spectrum Canon EOS 10D

Load measured solar spectra from files

i <- readFITS("Solar Spectrum Canon EOS 10D I.fit")
r <- readFITS("Solar Spectrum Canon EOS 10D R.fit")
g <- readFITS("Solar Spectrum Canon EOS 10D G.fit")
b <- readFITS("Solar Spectrum Canon EOS 10D B.fit")

Calculate spectra intensity

profile <- function(fits, band, scale) {
    means <- rowMeans(fits$imDat)
    data.frame(x    = seq(1:fits$axDat$len[1]), # X position in the spectrum
               y    = means / max(means),       # Normalized intensity
               raw  = means,                    # Raw intensity
               band = band                      # Colour band
    ) 
}

I <- profile(i, "I")
R <- profile(r, "R")
G <- profile(g, "G")
B <- profile(b, "B")

Adjust green values as we have two green units per the Bayer matrix element

G$raw <- 2 * G$raw

Renormalize red, green, and blue spectra

f <- max(R$raw, G$raw, B$raw)
R$y <- R$raw / f
G$y <- G$raw / f
B$y <- B$raw / f

Combine all spectra into one data frame

S <- data.frame(x=I$x, y=I$y, raw=I$raw, band=I$band)
for (df in list(R,G,B)) {
    S <- rbind(S, data.frame(x=df$x, y=df$y, raw=df$raw, band=df$band))
}

For wavelength calibration I have identified some lines in the spectra

calibration.data <- data.frame(
    x       = c(404,640,1072,1354),         # X position in the spectrum
    name    = c("G","H β", "D1+D2", "H α"), # Name of (Fraunhofer) line
    w       = c(430.7,486.1,589.3,653.3)    # wavelength [nm]
)
print(data.table(calibration.data))
##       x  name     w
## 1:  404     G 430.7
## 2:  640   H β 486.1
## 3: 1072 D1+D2 589.3
## 4: 1354   H α 653.3

Use a polynomial model to calibrate wavelength and calculate limits for plots

model <- lm(w ~ poly(x,3), calibration.data)
S$w <- predict(model, newdata=data.frame(x=S$x))

w.min <- max(min(S$w), 350)
w.max <- min(max(S$w), 750)
ggplot() +
    geom_line(data=subset(S,band=="I"), aes(x=w,y=y)) +
    geom_vline(data=calibration.data, aes(xintercept=w, color=c("1", "2", "3", "4"))) +
    scale_colour_manual(guide=FALSE, values=c("darkblue", "blue", "orange", "red")) +
    geom_text(data=calibration.data, aes(x=w+5,y=0.95,label=name), hjust=0) +
    theme_bw() +
    labs(x="Wavelength [nm]", y="Intensity", title="Solar Spectrum / Canon EOS 10D (I)") +
    coord_cartesian(xlim=c(w.min,w.max))

Load solar reference spectrum (AM0) from file

solar.reference <- read.csv("Wehrli85.csv", sep="", strip.white = T)
ggplot() +
    geom_line(data=solar.reference, aes(x=w,y=irradiance)) +
    theme_bw() +
    labs(x="Wavelength [nm]", y="Irradiance [W/m^2/nm]", title="Solar Reference Spectrum (Wehrli 1985)") +
    coord_cartesian(xlim=c(w.min,w.max))

Load atmosphere transmission from file

atmosphere <- read.csv("Atmosphere Transmission.csv", sep="")
ggplot() +
    geom_line(data=atmosphere, aes(x=w,y=100*transmission)) +
    theme_bw() +
    labs(x="Wavelength [nm]", y="Transmission [%]", title="Atmosphere Transmission (Airmass=2.13)") +
    coord_cartesian(xlim=c(w.min,w.max))

Calculate spectral response of Canon EOS 10D

S$irradiance <- approx(solar.reference$w, solar.reference$irradiance, S$w)$y
S$transmission <- approx(atmosphere$w, atmosphere$transmission, S$w)$y
S$response <- S$y / (S$irradiance * S$transmission)
smooth <- function(x, y, band, spar) {
    data.frame(x=x, y=smooth.spline(x, y, spar=spar)$y, band=band)
}

I <- subset(S, band=="I")
I$response <- I$response / max(I$response)

df <- smooth(I$w, I$response, I$band, 0.5)

ggplot() +
    geom_line(data=df, aes(x=x,y=y)) +
    theme_bw() +
    labs(x="Wavelength [nm]", y="Spectral Response", title="Spectral Response Canon EOS 10D (I)") +
    coord_cartesian(xlim=c(w.min,w.max))

R <- subset(S, band=="R")
G <- subset(S, band=="G")
B <- subset(S, band=="B")

f <- max(R$response, G$response, B$response)
R$response <- R$response / f
G$response <- G$response / f
B$response <- B$response / f

df <- smooth(R$w, R$response, R$band, 0.5)
df <- rbind(df, smooth(G$w, G$response, G$band, 0.5))
df <- rbind(df, smooth(B$w, B$response, B$band, 0.5))

ggplot() +
    geom_line(data=df, aes(x=x,y=y,color=band)) +
    scale_colour_manual(guide=FALSE, values=c("red", "green", "blue"), labels=c("R", "G", "B")) +
    theme_bw() +
    labs(x="Wavelength [nm]", y="Spectral Response", title="Spectral Reponse Canon EOS 10D (RGB)") +
    coord_cartesian(xlim=c(w.min,w.max))

The source is available as GitHub project.