//
//  Bifurcation.cpp
//  BasicVascularModeling
//
//  Created by Artur Klepaczko on 28.05.2014.
//  Copyright (c) 2014 Artur Klepaczko. All rights reserved.
//

#include "Bifurcation.h"
#include "Utils.h"
#include <math.h>

Bifurcation::Bifurcation(Vessel* bifurcatingVessel, Point3D newTerminalPoint) {
    
    double f1 = bifurcatingVessel->getFlow();
    double f2 = TerminalFlow;
    double f0 = f1 + f2;
    
//    printf("Bifurcation::Bifurcation(...). f1=%g, f2=%g, f0=%g.\n", f1, f2, f0);
    double x0 = bifurcatingVessel->getInlet().x;
    double x1 = bifurcatingVessel->getOutlet().x;
    double x2 = newTerminalPoint.x;
    double x_init = (f0*x0 + f1*x1 + f2*x2) / (2*f0);
    
    double y0 = bifurcatingVessel->getInlet().y;
    double y1 = bifurcatingVessel->getOutlet().y;
    double y2 = newTerminalPoint.y;
    double y_init = (f0*y0 + f1*y1 + f2*y2) / (2*f0);
    
    double z0 = bifurcatingVessel->getInlet().z;
    double z1 = bifurcatingVessel->getOutlet().z;
    double z2 = newTerminalPoint.z;
    double z_init = (f0*z0 + f1*z1 + f2*z2) / (2*f0);
    
    Point3D temporalBifurcationPoint;
    setPoint(&temporalBifurcationPoint, x_init, y_init, z_init);
    
    double l1 = distanceP2P(bifurcatingVessel->getOutlet(), temporalBifurcationPoint);
    
    double bifurcationPressure = bifurcatingVessel->getPressureOut() + f1*8*BloodViscosityInPa_s*l1/(M_PI*pow(bifurcatingVessel->getRadius(),4));
    
    if (bifurcatingVessel->hasParent()) {
        _parent = new Vessel(bifurcatingVessel->getInlet(), temporalBifurcationPoint, f0, bifurcatingVessel->getParent(), bifurcationPressure);
    }
    else {
        _parent = new Vessel(bifurcatingVessel->getInlet(), temporalBifurcationPoint, f0, bifurcatingVessel->getPressureIn(), bifurcationPressure);
    }
    
    _son = new Vessel(temporalBifurcationPoint, bifurcatingVessel->getOutlet(), f1, _parent, bifurcatingVessel->getPressureOut());
    _daughter = new Vessel(temporalBifurcationPoint, newTerminalPoint, f2, _parent, PressureAtOutlets);
    
    _parent->setChildren(_son, _daughter);
    
    _son->setRadius(bifurcatingVessel->getRadius());
    
    double l2 = distanceP2P(_daughter->getInlet(), _daughter->getOutlet());
    double r2 = pow(8*f2*BloodViscosityInPa_s*l2/(M_PI*(bifurcationPressure-PressureAtOutlets)), 0.25);
    _daughter->setRadius(r2);
    
    double r0 = radiusFromBifurcationLaw(_parent, _son, _daughter, BifurcationLawPower);
    _parent->setRadius(r0);
    
    _initVolume = bifurcatingVessel->getVolume();
}

Bifurcation::~Bifurcation() {
    _parent->DeleteVesel(true); // deletes all children automatically
}

double Bifurcation::bifurcationVolume() {
    
    double result = _parent->getVolume() + _son->getVolume() + _daughter->getVolume();
    return result;
}

void Bifurcation::optimize() {
    
    double current_volume = bifurcationVolume();
    double volumeIncrease = current_volume - _initVolume;
    double minIncrease = volumeIncrease;
    Point3D bestPoint = _parent->getOutlet();
    
    int N=10;

    double x0 = _parent->getInlet().x;
    double x1 = _son->getOutlet().x;
    double x2 = _daughter->getOutlet().x;
    
    double y0 = _parent->getInlet().y;
    double y1 = _son->getOutlet().y;
    double y2 = _daughter->getOutlet().y;
    
    double z0 = _parent->getInlet().z;
    double z1 = _son->getOutlet().z;
    double z2 = _daughter->getOutlet().z;

    double r0 = _parent->getRadius();
    double r1 = _son->getRadius();
    double r2 = _daughter->getRadius();
    
    double l0 = distanceP2P(_parent->getInlet(), _parent->getOutlet());
    double l1 = distanceP2P(_son->getInlet(), _son->getOutlet());
    double l2 = distanceP2P(_daughter->getInlet(), _daughter->getOutlet());
    
    double f1 = _son->getFlow();
    double f2 = _daughter->getFlow();
    
    for (int i=0; i<N; ++i) {
    
        double denominator = r0*r0/l0 + r1*r1/l1 + r2*r2/l2;
        
        double x_new = (x0*r0*r0/l0 + x1*r1*r1/l1 + x2*r2*r2/l2)/denominator;
        
        double y_new = (y0*r0*r0/l0 + y1*r1*r1/l1 + y2*r2*r2/l2)/denominator;
        
        double z_new = (z0*r0*r0/l0 + z1*r1*r1/l1 + z2*r2*r2/l2)/denominator;
        
        Point3D newTerminalPoint;
        setPoint(&newTerminalPoint, x_new, y_new, z_new);

        _parent->setOutllet(newTerminalPoint);
        _son->setInlet(newTerminalPoint);
        _daughter->setInlet(newTerminalPoint);
        
        l0 = distanceP2P(_parent->getInlet(), _parent->getOutlet());
        l1 = distanceP2P(_son->getInlet(), _son->getOutlet());
        l2 = distanceP2P(_daughter->getInlet(), _daughter->getOutlet());
        
        double bifurcationPressure = _son->getPressureOut() + f1*8*BloodViscosityInPa_s*l1/(M_PI*pow(r1,4));
        _parent->setPressureOut(bifurcationPressure);
        _son->setPressureIn(bifurcationPressure);
        _daughter->setPressureIn(bifurcationPressure);
        
        r2 = pow(8*f2*BloodViscosityInPa_s*l2/(M_PI*(bifurcationPressure-PressureAtOutlets)), 0.25);
        _daughter->setRadius(r2);
        
        r0 = radiusFromBifurcationLaw(_parent, _son, _daughter, BifurcationLawPower);
        _parent->setRadius(r0);
        
        current_volume = bifurcationVolume();
        volumeIncrease = current_volume - _initVolume;
        
//        printf("volume is %g (increase = %g).\n", current_volume, volumeIncrease);
        if (volumeIncrease < minIncrease) {
            minIncrease = volumeIncrease;
            bestPoint = _parent->getOutlet();
        }
    }
    
    _parent->setOutllet(bestPoint);
    _son->setInlet(bestPoint);
    _daughter->setInlet(bestPoint);
    
    l0 = distanceP2P(_parent->getInlet(), _parent->getOutlet());
    l1 = distanceP2P(_son->getInlet(), _son->getOutlet());
    l2 = distanceP2P(_daughter->getInlet(), _daughter->getOutlet());
    
    double bifurcationPressure = _son->getPressureOut() + f1*8*BloodViscosityInPa_s*l1/(M_PI*pow(r1,4));
    _parent->setPressureOut(bifurcationPressure);
    _son->setPressureIn(bifurcationPressure);
    _daughter->setPressureIn(bifurcationPressure);
    
    r2 = pow(8*f2*BloodViscosityInPa_s*l2/(M_PI*(bifurcationPressure-PressureAtOutlets)), 0.25);
    _daughter->setRadius(r2);
    
    r0 = radiusFromBifurcationLaw(_parent, _son, _daughter, BifurcationLawPower);
    _parent->setRadius(r0);
    
}