For what it is worth, I can see that the agreement with SNANA gets worse if we change the 
the linear interpolation of the error components to quadratic order
![snana_33c_sncosmo_modelerr_diag](https://cloud.githubusercontent.com/assets/2281358/2956953/67e449b2-da98-11e3-8b10-844ef53e9ac5.png)
Introducing the flexibility and doing M0, M1, Cl quadratically and the variance components linearly brings back the agreement. It is perhaps useful to have this flexibility.
![snana_33c_sncosmo_modelerr_diag](https://cloud.githubusercontent.com/assets/2281358/2957184/7feb68c6-da9b-11e3-99b3-b61d52a718a6.png)