Cops: C-OPS Processing

popt.d.e <- matrix(c(
	350.0,.778  ,355.0,.767  ,
	360.0,.756  ,365.0,.737  ,370.0,.720  ,375.0,.700  ,
	380.0,.685  ,385.0,.670  ,390.0,.655  ,395.0,.645  ,
	400.0,.64358,405.0,.64776,410.0,.65175,415.0,.65555,
	420.0,.65917,425.0,.66259,430.0,.66583,435.0,.66889,
	440.0,.67175,445.0,.67443,450.0,.67692,455.0,.67923,
	460.0,.68134,465.0,.68327,470.0,.68501,475.0,.68657,
	480.0,.68794,485.0,.68903,490.0,.68955,495.0,.68947,
	500.0,.68880,505.0,.68753,510.0,.68567,515.0,.68320,
	520.0,.68015,525.0,.67649,530.0,.67224,535.0,.66739,
	540.0,.66195,545.0,.65591,550.0,.64927,555.0,.64204,
	560.0,.6400,565.0,.6300,570.0,.6230,575.0,.6150,
	580.0,.6100,585.0,.6140,590.0,.6180,595.0,.6220,
	600.0,.6260,605.0,.6300,610.0,.6340,615.0,.6380,
	620.0,.6420,625.0,.6470,630.0,.6530,635.0,.6580,
	640.0,.6630,645.0,.6670,650.0,.6720,655.0,.6770,
	660.0,.6820,665.0,.6870,670.0,.6950,675.0,.6970,
	680.0,.6930,685.0,.6650,690.0,.6400,695.0,.6200,
	700.0,.6000
	), ncol = 2, byrow = T)
popt.f.e <- approxfun(popt.d.e, rule = 1)

popt.d.chi <- matrix(c(
	350.0,.153  ,355.0,.149  ,
	360.0,.144  ,365.0,.140  ,370.0,.136  ,375.0,.131  ,
	380.0,.127  ,385.0,.123  ,390.0,.119  ,395.0,.118  ,
	400.0,.11748,405.0,.12066,410.0,.12259,415.0,.12326,
	420.0,.12269,425.0,.12086,430.0,.11779,435.0,.11372,
	440.0,.10963,445.0,.10560,450.0,.10165,455.0,.09776,
	460.0,.09393,465.0,.09018,470.0,.08649,475.0,.08287,
	480.0,.07932,485.0,.07584,490.0,.07242,495.0,.06907,
	500.0,.06579,505.0,.06257,510.0,.05943,515.0,.05635,
	520.0,.05341,525.0,.05072,530.0,.04829,535.0,.04611,
	540.0,.04419,545.0,.04253,550.0,.04111,555.0,.03996,
	560.0,.0390,565.0,.0375,570.0,.0360,575.0,.0340,
	580.0,.0330,585.0,.0328,590.0,.0325,595.0,.0330,
	600.0,.0340,605.0,.0350,610.0,.0360,615.0,.0375,
	620.0,.0385,625.0,.0400,630.0,.0420,635.0,.0430,
	640.0,.0440,645.0,.0445,650.0,.0450,655.0,.0460,
	660.0,.0475,665.0,.0490,670.0,.0515,675.0,.0520,
	680.0,.0505,685.0,.0440,690.0,.0390,695.0,.0340,
	700.0,.0300
	), ncol = 2, byrow = T)
popt.f.chi <- approxfun(popt.d.chi, rule = 1)

# Pope et Fry a partir de 385
# de 350 a 375 Segandares et Fry
# 380 n'est la valeur de Pope et Fry : c'est une valeur de raccordement
popt.d.aw <- matrix(c(
	350.0,0.02037, 355.0,0.01819, 360.0,0.01565, 365.0,0.01319,
	370.0,0.01241, 375.0,0.01096,
	380.0,0.01040, 382.5,0.01044, 385.0,0.00941, 387.5,0.00917,
	390.0,0.00851, 392.5,0.00829,
	395.0,0.00813, 397.5,0.00775, 400.0,0.00663, 402.5,0.00579,
	405.0,0.00530, 407.5,0.00503,
	410.0,0.00473, 412.5,0.00452, 415.0,0.00444, 417.5,0.00442,
	420.0,0.00454, 422.5,0.00474,
	425.0,0.00478, 427.5,0.00482, 430.0,0.00495, 432.5,0.00504,
	435.0,0.00530, 437.5,0.00580,
	440.0,0.00635, 442.5,0.00696, 445.0,0.00751, 447.5,0.00830,
	450.0,0.00922, 452.5,0.00969,
	455.0,0.00962, 457.5,0.00957, 460.0,0.00979, 462.5,0.01005,
	465.0,0.01011, 467.5,0.0102,
	470.0,0.0106, 472.5,0.0109, 475.0,0.0114, 477.5,0.0121,
	480.0,0.0127, 482.5,0.0131,
	485.0,0.0136, 487.5,0.0144, 490.0,0.0150, 492.5,0.0162,
	495.0,0.0173, 497.5,0.0191,
	500.0,0.0204, 502.5,0.0228, 505.0,0.0256, 507.5,0.0280,
	510.0,0.0325, 512.5,0.0372,
	515.0,0.0396, 517.5,0.0399, 520.0,0.0409, 522.5,0.0416,
	525.0,0.0417, 527.5,0.0428,
	530.0,0.0434, 532.5,0.0447, 535.0,0.0452, 537.5,0.0466,
	540.0,0.0474, 542.5,0.0489,
	545.0,0.0511, 547.5,0.0537, 550.0,0.0565, 552.5,0.0593,
	555.0,0.0596, 557.5,0.0606,
	560.0,0.0619, 562.5,0.0640, 565.0,0.0642, 567.5,0.0672,
	570.0,0.0695, 572.5,0.0733,
	575.0,0.0772, 577.5,0.0836, 580.0,0.0896, 582.5,0.0989,
	585.0,0.1100, 587.5,0.1220,
	590.0,0.1351, 592.5,0.1516, 595.0,0.1672, 597.5,0.1925,
	600.0,0.2224, 602.5,0.2470,
	605.0,0.2577, 607.5,0.2629, 610.0,0.2644, 612.5,0.2665,
	615.0,0.2678, 617.5,0.2707,
	620.0,0.2755, 622.5,0.2810, 625.0,0.2834, 627.5,0.2904,
	630.0,0.2916, 632.5,0.2995,
	635.0,0.3012, 637.5,0.3077, 640.0,0.3108, 642.5,0.322,
	645.0,0.325, 647.5,0.335,
	650.0,0.340, 652.5,0.358, 655.0,0.371, 657.5,0.393,
	660.0,0.410, 662.5,0.424,
	665.0,0.429, 667.5,0.436, 670.0,0.439, 672.5,0.448,
	675.0,0.448, 677.5,0.461,
	680.0,0.465, 682.5,0.478, 685.0,0.486, 687.5,0.502,
	690.0,0.516, 692.5,0.538,
	695.0,0.559, 697.5,0.592, 700.0,0.624, 702.5,0.663,
	705.0,0.704, 707.5,0.756,
	710.0,0.827, 712.5,0.914, 715.0,1.007, 717.5,1.119,
	720.0,1.231, 722.5,1.356,
	725.0,1.489, 727.5,1.678, 775.0,2.400, 865.0,5.550
	), ncol = 2, byrow = T)
popt.f.aw <- approxfun(popt.d.aw, rule = 1)

popt.d.mud <- list(
	chlor = c(0.03,0.1,0.3,1.,3.,10.),
	w = c(350,400,412,443,490,510,555,620,670,700),
	chlormin = 0.03,
	chlormax = 10,
	wmin = 350,
	wmax = 700,
	val = matrix(c(
		0.770,0.769,0.766,0.767,0.767,0.767,
		0.770,0.769,0.766,0.767,0.767,0.767,
		0.765,0.770,0.774,0.779,0.782,0.782,
		0.800,0.797,0.796,0.797,0.799,0.799,
		0.841,0.824,0.808,0.797,0.791,0.791,
		0.872,0.855,0.834,0.811,0.796,0.796,
		0.892,0.879,0.858,0.827,0.795,0.795,
		0.911,0.908,0.902,0.890,0.871,0.871,
		0.914,0.912,0.909,0.901,0.890,0.890,
		0.914,0.912,0.909,0.901,0.890,0.890
		),
		nrow = 10,
		ncol = 6,
		byrow = T)
)
dimnames(popt.d.mud$val) <- list(popt.d.mud$w, popt.d.mud$chlor)
popt.f.mud <- function(wl,chl) {
	if(wl < popt.d.mud$wmin) return(NA)
	if(wl > popt.d.mud$wmax) return(NA)
	if(chl < popt.d.mud$chlormin) return(NA)
	if(chl > popt.d.mud$chlormax) return(NA)
	i = max(which(popt.d.mud$w <= wl))
	i = min(i,9)
	j = max(which(popt.d.mud$chlor <= chl))
	j = min(j,5)
	mud1 <- approx(popt.d.mud$w[i:(i+1)],popt.d.mud$val[i:(i+1),j    ],wl)$y
	mud2 <- approx(popt.d.mud$w[i:(i+1)],popt.d.mud$val[i:(i+1),j + 1],wl)$y
	return(approx(popt.d.mud$chlor[j:(j+1)], c(mud1,mud2), chl)$y)
}

popt.f.a <- function(wl,chl) {
	chi=popt.f.chi(wl)
	e=popt.f.e(wl)
	mud=popt.f.mud(wl,chl)

	bbw=0.5*popt.f.bw(wl)
	if(abs(chl) < 0.0001) {
		bbp=0.
		kd=popt.f.kw(wl)
	} else {
		bp=(popt.f.b550(chl)-popt.f.bw(550.))*popt.f.sdobp(wl,chl)
		bbpt=0.002+0.01*(0.5-0.25*log10(chl))
		bbp=bbpt*bp
		kd=popt.f.kw(wl)+chi*chl^e
	}
	bb=bbw+bbp

	u1=0.75
	r1=0.33*bb/u1/kd
	err <- 1
	while(err >= 0.0001) {
		u2=mud*(1.-r1)/(1.+mud*r1/0.42)
		r2=0.33*bb/u2/kd
		err=abs((r2-r1)/r2)
		r1=r2
	}
	return(u2*kd)
}

popt.f.b550 <- function(chl) {
      return(0.416*chl^0.767)
}

popt.f.bw <- function(x) {
	return(0.00288/(x/500.)^4.32)
}

popt.f.kw <- function(x) {
      return(popt.f.aw(x)+0.5*popt.f.bw(x))
}

popt.f.ref <- function(wl,chl) {
	chi=popt.f.chi(wl)
	e=popt.f.e(wl)
	mud=popt.f.mud(wl,chl)

	bbw=0.5*popt.f.bw(wl)
	if(abs(chl) < 0.0001) {
		bbp=0.
		kd=popt.f.kw(wl)
	} else {
		bp=(popt.f.b550(chl)-popt.f.bw(550.))*popt.f.sdobp(wl,chl)
		bbpt=0.002+0.01*(0.5-0.25*log10(chl))
		bbp=bbpt*bp
		kd=popt.f.kw(wl)+chi*chl^e
	}
	bb=bbw+bbp

	u1=0.75
	r1=0.33*bb/u1/kd
	err <- 1
	while(err >= 0.0001) {
		u2=mud*(1.-r1)/(1.+mud*r1/0.42)
		r2=0.33*bb/u2/kd
		err=abs((r2-r1)/r2)
		r1=r2
	}
	return(r2)
}

popt.f.sdobp <- function(wl,chl) {
      if(chl < 2) {
        expo = -0.5 * log10(0.5 * chl)
      } else {
        expo = 0.
      }
      return((550./wl)^expo)
}

popt.f.f <- function(x0) {
	if(length(x0) != 3) stop("erreur calcul f: length(argument) != 3")
	#table.f <- list(w , ts , lchl , f )
	x <- numeric(0)
	dx <- array(NA, dim = c(2,3))
	ix <- integer(0)
	ip <- integer(0)
	for(i in 1:3) {
		x[i] <- min(max(x0[i],min(table.f[[i]])),max(table.f[[i]]))
		di <- table.f[[i]] - x[i]
		ind <- which(di > 0, arr.ind=TRUE)
		if(length(ind) > 0) ix[i] <- ind[1]
		else ix[i] <- which(di == 0, arr.ind=TRUE)
		ip[i] <- ix[i]-1
		dx[1,i] <- di[ix[i]]
		dx[2,i] <- -di[ip[i]]
		dx[,i] <- dx[,i] / sum(dx[,i])
	}
	tab <- table.f$f[ip[1]:ix[1],ip[2]:ix[2],ip[3]:ix[3]]
	y <- 0
	for (i in 1:2) {
		a = dx[i,1]
		for (j in 1:2) {
			b = a * dx[j,2]
			for (k in 1:2) {
				y = y + b * dx[k,3] * tab[i,j,k]
			}
		}
	}
	return( list(w = x[1], ts = x[2], lchl = x[3], f =y))
}

popt.f.Q <- function(x0) {
	if(length(x0) != 5) stop("erreur calcul Q: length(argument) != 5")
	#table.Q <- list(w , ts , lchl , tv , phi , Q )
	x <- numeric(0)
	dx <- array(NA, dim = c(2,5))
	ix <- integer(0)
	ip <- integer(0)
	for(i in 1:5) {
		x[i] <- min(max(x0[i],min(table.Q[[i]])),max(table.Q[[i]]))
		di <- table.Q[[i]] - x[i]
		ind <- which(di > 0, arr.ind=TRUE)
		if(length(ind) > 0) ix[i] <- ind[1]
		else ix[i] <- which(di == 0, arr.ind=TRUE)
		ip[i] <- ix[i]-1
		dx[1,i] <- di[ix[i]]
		dx[2,i] <- -di[ip[i]]
		dx[,i] <- dx[,i] / sum(dx[,i])
	}
	tab <- table.Q$Q[ip[1]:ix[1],ip[2]:ix[2],ip[3]:ix[3],ip[4]:ix[4],ip[5]:ix[5]]
	y <- 0
	for (i in 1:2) {
		a = dx[i,1]
		for (j in 1:2) {
			b = a * dx[j,2]
			for (k in 1:2) {
				c = b * dx[k,3]
				for (l in 1:2) {
					d = c * dx[l,4]
					for (m in 1:2) {
						y = y + d * dx[m,5] * tab[i,j,k,l,m]
					}
				}
			}
		}
	}
	return( list(w = x[1], ts = x[2], lchl = x[3], tv = x[4], phi = x[5], Q =y))
}

################ fausse vectorisation ##############
popt.fv.ref <- function(wl,chl) {
	r <- array(NA, dim = c(length(wl),length(chl)), dimnames = list(wl,chl))
	i = 0
	for (w in wl) {
		i = i + 1
		j = 0
		for (c in chl) {
			j = j + 1
			r[i,j] <- popt.f.ref(w,c)
		}
	}
	return(r)
}
popt.fv.a <- function(wl,chl) {
	r <- array(NA, dim = c(length(wl),length(chl)), dimnames = list(wl,chl))
	i = 0
	for (w in wl) {
		i = i + 1
		j = 0
		for (c in chl) {
			j = j + 1
			r[i,j] <- popt.f.a(w,c)
		}
	}
	return(r)
}