// Federico Ziraldo 92909
// Daniel Madrid 95483

#include <iostream>

#define D -0.14
#define U 39 // km/h

using namespace std;

void resolucion(long double o, long double K, long double h, long double xn,
 long double cn, long double dcn, long double* tablaCn, long double* tablaXn){
    int i = 1;
    while (cn > .1){
        long double dcnmenos1 = dcn;
        xn += h;
        dcn = (dcn + h/D*(1-o)*(U*dcn + K*cn) + h*o*K/D*cn + h*h*o*K/D*(1-o)*dcn) / (1 - h*o*U/D - h*h*o*o*K/D);
        cn = cn + h*((1-o)*dcnmenos1 + o*dcn);
        cout << "error: " << (tablaCn[int(h*10)*i]-cn) / tablaCn[int(h*10)*i] << "    "
         << xn << "    " << cn << "    " << tablaCn[int(h*10)*i] << "    " << tablaXn[int(h*10)*i] << endl;
        ++i;
    }
}

int main()
{
    long double tablaXn[1100];
    long double tablaCn[1100];
    long double K = 10, o, h = .1, xn = 0,
     cn = 1000000, // mg/m3
     dcn = -1000; // mg/(m3*km)
    cout << "ingresar o. h = .1" << endl;
    cin >> o ;
    int i = 0;
    tablaXn[i] = xn;
    tablaCn[i] = cn;
    while (cn > .1 || xn < 101){
        ++i;
        long double dcnmenos1 = dcn;
        xn += h;
        dcn = (dcn + h/D*(1-o)*(U*dcn + K*cn) + h*o*K/D*cn + h*h*o*K/D*(1-o)*dcn) / (1 - h*o*U/D - h*h*o*o*K/D);
        cn = cn + h*((1-o)*dcnmenos1 + o*dcn);
        tablaXn[i] = xn;
        tablaCn[i] = cn;
    }

    xn = 0;
    cn = 1000000; // mg/m3
    dcn = -1000; // mg/(m3*km)
    cout << endl;
    resolucion(o, K, .4, xn, cn, dcn, tablaCn, tablaXn);

    return 0;
}