source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/G4fresnellens/src/ParabolicSurface.cc @ 114

#include "ParabolicSurface.h"
#include "LensUtils.h"

using std::swap;
using namespace G4FresnelLens;

//______________________________________________________________________________
ParabolicSurface::ParabolicSurface(double k, double z, double rhomax):
ParametricSurface(z,rhomax),
fK(k)
{
    //kHalfTolerance=kTolerance;
    //kTolerance=2.*kHalfTolerance;
    Configure();
}

//______________________________________________________________________________
ParabolicSurface::~ParabolicSurface(){


}

//______________________________________________________________________________
void ParabolicSurface::Configure(){
    fZ1=fZmin=PointOnSurf(0.);
    fZ2=fZmax=PointOnSurf(fR);
    if ( fZmin>fZmax ) swap(fZmin,fZmax);
}

//__________________________________________________________________________________
void ParabolicSurface::print(){
    printf("paraboloid z=%g perp=%g k=%g min/max=%g/%g\n", fZpos, fR, fK, fZmin, fZmax);
}

//__________________________________________________________________________________
void ParabolicSurface::GuessParameters(int n, int tn, int* idx, double* rho, double* z, bool minimum){
    if ( fK!=DBL_MAX ) return;
    int i1a=idx[1];
    int i1b=i1a+1;
    int i1=(z[i1a]>z[i1b])^minimum?i1a:i1b;

    int i2a=idx[tn-1];
    int i2b=i2a+1;
    int i2=(z[i2a]>z[i2b])^minimum?i2a:i2b;

    double dz=z[i2]-z[i1];
    double dr2=rho[i2]*rho[i2]-rho[i1]*rho[i1];
    fK=dr2?dz/dr2:0.;
}

//__________________________________________________________________________________
void ParabolicSurface::GuessParameters(const std::vector<SSurface*>& surfaces, bool minimum){
    if ( fK!=DBL_MAX ) return;

    SSurface* s1=*(surfaces.begin());
    SSurface* s2=*(surfaces.end()-1);
    double r1,r2,z1,z2;
    double za=s1->GetZ1();
    double zb=s1->GetZ2();
    if ( (za>zb)^minimum ){
        z1=s1->GetZ1();
        r1=s1->GetR1();
    }
    else {
        z1=s1->GetZ2();
        r1=s1->GetR2();
    }

    za=s2->GetZ1();
    zb=s2->GetZ2();
    if ( (za>zb)^minimum ){
        z2=s2->GetZ1();
        r2=s2->GetR1();
    }
    else {
        z2=s2->GetZ2();
        r2=s2->GetR2();
    }

    double dz=z2-z1;
    double dr2=r2*r2-r1*r1;
    fK=dr2?dz/dr2:0.;
}

//__________________________________________________________________________________
double ParabolicSurface::DistanceToSurf(const G4ThreeVector& p, const G4ThreeVector& v){
    double perp2=p.perp2();
    double vperp2=v.perp2();
    double a = fK*vperp2;
    double b = 2.*fK*(p.x()*v.x()+p.y()*v.y())-v.z();
    double c = fK*perp2-p.z()+fZpos;
    double d1, d2;
    int nroot=SolveQuadratic(a, b, c, d1, d2);

    if ( !nroot ) return kInfinity;
    #ifdef DEBUG
        printf("pos(%.10g, %.10g, %.10g), dir(%.10g, %.10g, %.10g)\n", p.x(), p.y(), p.z(), v.x(), v.y(), v.z());
        printf("abc: %.10g, %.10g, %.10g\n", a,b,c);
        double d=sqrt(b*b-4*a*c);
        double x1=(-b+d)/2./a;
        printf("s(d)=%.10g, d1:d2: %.10g, %.10g\n",d, (-b-d)/2./a, x1);
        printf("s(d)=%.10g, d1:d2: %.10g, %.10g\n",d, x1, -b/a-x1);
        printf("%i roots: %.10g, %.10g\n", nroot, d1, d2);
    #endif
    G4ThreeVector pos;
    if ( d1>=0. ) {
        pos = p+v*d1;
        double zpos=pos.z();
        double perpa2=pos.perp2();
        if ( perpa2<=fR*fR && zpos>=fZmin && zpos<=fZmax ) {
            return d1;
        }
    }
    if ( nroot==1 || d2<0. ) return kInfinity;

    pos = p+v*d2;
    double zpos=pos.z();
    double perpa2=pos.perp2();
    if ( perpa2<=fR*fR && zpos>=fZmin && zpos<=fZmax ) {
        return d2;
    }
    return kInfinity;
}

//__________________________________________________________________________________
double ParabolicSurface::DistanceToSurfSafe(const G4ThreeVector& p){
    //
    // return undirectional distance to surface
    // can be underestimated
    // thus returning distance to the bounding tube
    //
    double perp2=p.perp2();
    double z=p.z();

    if ( perp2<=fR*fR ) {
        if ( z>fZmax ) return z-fZmax;
        if ( z<fZmin ) return fZmin-z;
        return 0.;
    }
    double perp=sqrt(perp2);
    if ( z>fZmax ) return Hyp(perp, z, fR, fZmax);
    if ( z<fZmin ) return Hyp(perp, z, fR, fZmin);
    return perp-fR;
}