#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem Flat_EarthDrag
#        ODE of dimension 8
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers
#
#       Authors:
#
#       Jeff R. Cash 
#            (Department of Mathematics, Imperial College,  London, England.)
#       Davy  Hollevoet 
#            (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium.)
#       Francesca Mazzia  
#            (Dipartimento di Matematica, Universita' di Bari, Italy)
#       Abdelhameed Nagy Abdo
#            (Dipartimento di Matematica, Universita' di Bari, Italy)
#            (Dept. of Mathematics, Faculty of Sciences, Benha  University,Egypt)


Flat_EarthDrag<- function(){
  
  
prob<- function(){ 
          fullnm <- 'Problem Flat_EarthDrag'
          problm = 'Flat_EarthDrag'
          typebvp   = 'BVP'
          
          neqn = 8
          nlbc = 4
          aleft    = 0
          aright   = 1
          numjac   = TRUE
          numbcjac = TRUE
          linear   = FALSE
          Rpar<-c(0.5)
          Ipar<-0
          return(list(fullnm=fullnm,problm=problm,typebvp=typebvp,neqn=neqn,
                      nlbc=nlbc,aleft=aleft,aright=aright,Rpar=Rpar,Ipar=Ipar,
                      numjac=numjac,numbcjac=numbcjac,linear=linear))
}
  
  
  
init <- function(neqn) { 
        givmsh = TRUE
        givey = TRUE
        nmsh = 11
        
        yval <- rep(0,neqn)
        yval[6]=-1
        yval[neqn]=100
        
        xguess = seq(0,1,by=1/(nmsh-1))
        yguess <- NULL
        for (i in 1:neqn)
                   yguess <- rbind(yguess,rep(yval[i],nmsh))
        return(list(givmsh=givmsh,givey=givey,
                    nmsh=nmsh,xguess=xguess,yguess=yguess))
}
  
feval = function(tau,X,eps,Rpar,Ipar){ 
          f<-2.1*10^6
          h<-180000
          m<-60880
          g_accel<-9.80665
          vc<-1000*sqrt((3.986004*10^5)/(6378.14+(h/1000)))
          beta<-180000/840
          eta<-1.225*0.5*7.069/2
    
          xbardot <- X[3]*(vc/h)
          ybardot <- X[4]*(vc/h)
          Vxbardot <- (f/vc*(-X[6]/sqrt(X[6]^2+X[7]^2))-eta*exp(-X[2]*beta)*X[3]*sqrt(X[3]^2+X[4]^2)*vc)/m
          Vybardot <- (f/vc*(-X[7]/sqrt(X[6]^2+X[7]^2))-eta*exp(-X[2]*beta)*X[4]*sqrt(X[3]^2+X[4]^2)*vc)/m-g_accel/vc
    if (sqrt(X[3]^2 + X[4]^2) == 0){
    
                                    lambda_2_bar <- 0
                                    lambda_3_bar <- 0
                                    lambda_4_bar <- -X[5]*(vc/h) 
                                    } else {
                                    lambda_2_bar <- -(X[6]*X[3]+X[7]*X[4])*eta*beta*sqrt(X[3]^2+X[4]^2)*exp(-X[2]*beta)*vc/m
                                    lambda_3_bar <- eta*exp(-X[2]*beta)*vc*(X[6]*(2*X[3]^2+X[4]^2)+X[7]*X[3]*X[4])/sqrt(X[3]^2+X[4]^2)/m
                                    lambda_4_bar <- -X[5]*vc/h+eta*exp(-X[2]*beta)*vc*(X[7]*(X[3]^2)+2*X[4]^2+X[6]*X[3]*X[4])/sqrt(X[3]^2+X[4]^2)/m}
    tf <- X[8];
    
    
    
            f<-(tf*c(xbardot, 
                      ybardot, 
                      Vxbardot,
                      Vybardot,
                      lambda_2_bar,
                      lambda_3_bar,
                      lambda_4_bar,
                      0))
  

            return(list(f))
}
  


jeval = function(x,y,eps,Rpar,Ipar){
        dfy <- matrix( nrow=6,ncol=6,byrow = TRUE, data = c(
         ))
    
        return((dfy))             
}
  
  
bcfun <- function(i, y, eps,Rpar,Ipar) {
          f<-2.1*10^6
          
          h<-180000
          m<-60880
          g_accel<-9.80665
          vc<-1000*sqrt((3.986004*10^5)/(6378.14+(h/1000)))
          beta<-180000/840
          eta<-1.225*0.5*7.069/2
          
          if (i == 1) return(y[1])
          if (i == 2) return(y[2])
          if (i == 3) return(y[3])
          if (i == 4) return(y[4])
          if (i == 5) return(y[2]-1.0e0)
          if (i == 6) return(y[3]-1.0e0)
          if (i == 7) return(y[4])
          if (i == 8) return ((-sqrt(y[6]^2+y[7]^2)*f/m/vc-(y[6]*y[3])*eta*exp(-beta)*sqrt(y[3]^2)*vc/m-y[7]*g_accel/vc)*y[8]+ 1)
    
  }
  
  dbcfun <- function(i, y, eps,Rpar,Ipar) {
    }
  
  
  
setoutput<- function(neqn,plotsol=NULL){ 
    solref = FALSE
    
    if (is.null(plotsol)){ 
                           nindsol = neqn
                           indsol = 1;neqn
                          }
    else{
                           nindsol = length(plotsol)
                           indsol = plotsol
                          }         
    return(list(solref=solref,nindsol=nindsol,indsol=indsol))
}
  
esolu <- function(X,parms,Rpar,Ipar){
    lambda=parms
    Exact =matrix( c( exp(-X/(lambda)),
                      (-exp(-X/(lambda)))/(lambda)),
                      ncol=2,nrow=length(X),byrow=FALSE)
    return(Exact=Exact)

    }
  
  return(list(prob=prob,init=init,feval=feval,jeval=jeval,bcfun=bcfun,dbcfun=dbcfun,setoutput=setoutput,esolu=esolu))   
  
}