anchor.diff:          if area = 0.0 then
               None
           else
  -            //max area is considered to be ~ 20
  -            //should give velocity of 1.0
  -            //thus
  -            let k = v_max * (1.0 - v_p) / max_hole
  -            let b = v_max * v_p
  -            let velocity = k * area + b                    
  -            #if DEBUG
  -            assert(abs(k* max_hole + b - v_max)<1e-14)
  -            assert(abs(k*0.0 + b - v_max*v_p)<1e-14)
  -            #endif
  -            let k = d_max * (1.0 - d_p) / max_hole
  -            let b = d_max * d_p
  -            let density = k * area + b                    
  -            #if DEBUG
  -            assert(abs(k* max_hole + b - d_max)<1e-14)
  -            assert(abs(k*0.0 + b - d_max*d_p)<1e-14)
  -            #endif            
  +            let v2 = v3p*v2p
  +            let d2 = d3p*d2p
  +            let va,vb,vc = get2ndPolyCoefs 0.0 (max_hole*vap) max_hole (v2*v1p) v2 v3p
  +            let da,db,dc = get2ndPolyCoefs 0.0 (max_hole*dap) max_hole (d2*d1p) d2 d3p
  +            let areaSq = area*area
               let burst = {
                   EmergenceTime=emergence_tick;
  -                Velocity=velocity;
  -                Density= density
  +                Velocity= va*areaSq + vb*area + vc;
  +                Density= da*areaSq + db*area + dc
               }
               Some(burst)
       Table.Map ["ts";"Relative_CH_CorrectSphereArea_j"] generatePulse predictorsT |> Seq.choose (fun b -> b)  |> List.ofSeq
  
positive.code: //             |> Array.Parallel.map current_lglk
               |> Array.map current_lglk
               |> Array.sum         
      if res>bestLglk then
          printfn "\nlglk improvement:%g\tVm:%g\tVp:%g\tDm:%g\tDp:%g\tVbg:%g\tDbg:%g\tSig:%g\tD:%g"
              res v_max v_p d_max d_p v_b d_b Sigma D
          bestLglk <- res
      else
          printf "."
      res
          
  
  [<EntryPoint>]
  let main argv =     
      //config
      let doEstimate = true
      let thinnObs = false
  
negative.code: 
      let estimate predictorsT observationsT =
          let parameters = Parameters.Empty
                              .Add("V_p3",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("V_p2",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("V_p1",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("D_p3",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("D_p2",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("D_p1",[|0.8|],eps,1.0,isLog=true,delay=0)
                              .Add("D_p",[|0.5|],eps,1.0,isLog=true,delay=0)
                              .Add("Area_vp",[|0.5|],eps,1.0,isLog=true,delay=0)
                              .Add("Area_dp",[|0.5|],eps,1.0,isLog=true,delay=0)
                              .Add("V_bg",[|0.1|],eps,0.5,isLog=true,delay=0)
                              .Add("D_bg",[|0.5|],eps,1.0,isLog=true,delay=0)
                              //.Add("D_bg",[|0.5|])
                              .Add("D",[|1500.0|],1200.0,1800.0,isLog=false,delay=0)
                              .Add("Sigma",[|0.01|],eps, 0.5,isLog=true,delay=0)                            
  