Changeset 19597
 Timestamp:
 Sep 2, 2019, 8:12:27 PM (6 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Mars/msimreflector/MFresnelLens.cc
r19595 r19597 193 193 const double c = n*u.fVectorPart; 194 194 195 const double R = 1  r*r*(1c*c); 196 197 const auto v = r*c + sqrt(R<0 ? 0 : R); 195 const double rc = r*c; 196 197 const double R = 1  r*r + rc*rc; // = 1  r*r *(1c*c); 198 199 const double v = rc + sqrt(R<0 ? 0 : R); 198 200 199 201 u = r*u  n*v; … … 202 204 Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const 203 205 { 204 const double F = 50; // [mm] focal length 205 const double R = 25; // [mm] radius of lens 206 const int N = 25; // [cnt] Nuber of grooves within radius 207 208 const double w = R/N; // [mm] Width of a single groove 209 210 const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength); 206 // Corsika Coordinates are in cm! 207 208 const double D = 54.92; // [cm] Diameter 209 const double F = 50.21; // [cm] focal length (see also MGeomCamFAMOUS!) 210 211 const double R = D/2; // [cm] radius of lens 212 const double w = 0.0100; // [cm] Width of a single groove 213 const double H = 0.18; // [cm] Thickness 214 215 // Focal length is given at this wavelength 216 const double n0 = CalcRefractiveIndex(546); 217 218 const double n = CalcRefractiveIndex(wavelength==0 ? 546 : wavelength); 219 220 // Ray has missed the lens 221 if (p.R()>R) 222 return 1; 223 224 const double phi = atan(p.R()/F); 211 225 212 226 const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove 213 227 214 const double r0 = nth*w;// Lower radius of nth groove215 const double r1 = (nth+1)*w;// Upper radius of nth groove216 const double cen= (nth+0.5)*w; // Center of nth groove228 const double r0 = nth *w; // Lower radius of nth groove 229 const double r1 = (nth+1) *w; // Upper radius of nth groove 230 const double rc = (nth+0.5)*w; // Center of nth groove 217 231 218 232 // FIXME: Do we have to check the nth1 and nth+1 groove as well? 219 233 220 // Angle to the center of the groove 221 const double phi = acos(cen / (F/1.125)); 222 223 // theta peak_z 224 const Groove g( TMath::Pi()/2  phi, r0*tan(phi)); 234 // Angle of grooves 235 // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf 236 237 // This might be a better reference (Figure 1 Variant 1) 238 // https://link.springer.com/article/10.3103/S0003701X13010064 239 240 // I suggest we just do the calculation ourselves (but we need to 241 // find out how our lens looks like) 242 // Not 100% sure, but I think we have SC943 243 // https://www.orafol.com/en/europe/products/opticsolutions/productlines#pl1 244 245 // sin(omega) = R / sqrt(R^2+f^2) 246 // tan(alpha) = sin(omega) / [ 1  sqrt(n^2sin(omega)^2) ] 247 248 const double so = rc / sqrt(rc*rc + F*F); 249 const double theta = atan(so / (1sqrt(n0*n0  so*so))); 250 251 // theta peak_z 252 const Groove g(TMath::Pi()/2  theta, r0*tan(theta)); 225 253 226 254 // Calculate the insersection between the ray and the cone 227 float dz = CalcIntersection(p, u, g,  3);255 float dz = CalcIntersection(p, u, g, H); 228 256 229 257 // No groove was hit … … 244 272 // FIXME: Surface roughness? 245 273 TVector3 norm; 246 norm.SetMagThetaPhi(1, g.fTheta +TMath::Pi()/2, p.XYvector().Phi());274 norm.SetMagThetaPhi(1, g.fThetaTMath::Pi()/2, p.XYvector().Phi()); 247 275 248 276 // Apply refraction at lens entrance (change directional vector) … … 253 281 const double v = u.fRealPart; // c changes in the medium 254 282 u.fRealPart = n*v; 255 p.PropagateZ(u,  3);283 p.PropagateZ(u, H); 256 284 u.fRealPart = v; // Set back to c 257 285 258 // New normal surface (bottom of lens)259 norm.SetMagThetaPhi(1, 0, 0);260 261 286 // Apply refraction at lens exit (change directional vector) 287 // Normal surface (bottom of lens) 262 288 // FIXME: Surace roughness 263 ApplyRefraction(u, norm, n, 1); 264 265 // To make this consistent with a mirror system, we now change our coordinate system 289 ApplyRefraction(u, TVector3(0, 0, 1), n, 1); 290 291 // To make this consistent with a mirror system, 292 // we now change our coordinate system 266 293 // Rays from the lens to the camera are upgoing (positive sign) 267 294 u.fVectorPart.SetZ(u.Z()); 268 295 269 // This is only correct if the focal distance is calculated from the inner lens surface! 270 p.fVectorPart.SetZ(0); 296 // In the datasheet, it looks as if F is calculated 297 // towards the center of the lens. 298 p.fVectorPart.SetZ(H/2); 271 299 272 300 return nth; // Keep track of groove index
Note: See TracChangeset
for help on using the changeset viewer.