Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Introduction

I created the following code to simulate the many possible interactions between a cylindrical protein crystal and a x-ray beam during serial femtosecond crystallography.

Manipulate[
  Graphics3D[
    {
     { EdgeForm[None], Directive[Yellow, Opacity[0.4], Specularity[White, 20]], 
       GeometricTransformation[Cylinder[{{-h + x, y, z}, {h + x, y, z}}, 1*10^-6], 
       RotationTransform[α Pi/2, {1, 0, 0}].RotationTransform[β Pi/2, {0, 1, 0}]]},
     { EdgeForm[None], Directive[Blue, Opacity[0.3], Specularity[White, 20]], 
       Cylinder[{{-1*10^-5, 0, 0}, {1*10^-5, 0, 0}}, r]} },

    Axes -> True, AxesLabel -> {"X", "Y", "Z"} ],

  {{α, 0}, -1, 1}, {{β, 1}, -1, 1}, 
  {{r, 1*10^-7}, 1*10^-7, 6*10^-6} , 
  {{h, 5*10^-6}, 1*10^-6, 10*10^-6},
  {{x, 0}, -1*10^-5, 1*10^-5},
  {{y, 0}, -1*10^-5, 1*10^-5}, 
  {{z, 0}, -1*10^-5, 1*10^-5} ]

Problem

I would like to be able to display the volume of the intersection of the two objects. As far as I have looked, this doesn't seem to be possible using only Cylinder and Graphics3D. If there is a way to find this volume using only Graphics objects that someone knows of, please enlighten me.

Also, the use of RegionPlot3D and NIntegrate as demonstrated by this example looks promising.

RegionPlot3D[x^2 + y^2 < 1, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]
NIntegrate[Boole[x^2 + y^2 < 1 ], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]

But my problem with using these functions is I don't know how to apply the Manipulate function to these raw inequalities in order to get the geometric transforms, variable height, translation, and change in beam radius like I have in my Graphics3D model.

Really all I care about it that it is manipulatable like the Graphics3D model I have shown and that the volume of the intersection is displayed as the model is manipulated. I greatly appreciate any help on this.

Solution

cyl[{x_, y_, z_}, r_, h_] := y^2 + z^2 < r^2 && -h < x < h
xform[x_, y_, z_, a_, b_] := RotationTransform[a Pi/2, {1, 0, 0}].RotationTransform[b Pi/2, {0, 1, 0}].TranslationTransform[{x, y, z}]

Manipulate[
 Graphics3D[
  {
   {EdgeForm[None], Directive[Yellow, Opacity[0.4], Specularity[White, 20]],
   GeometricTransformation[ Cylinder[{{-h, 0, 0}, {h, 0, 0}}, 1*10^-6], xform[x, y, z, a, b]]},
   {EdgeForm[None], Directive[Blue, Opacity[0.3], Specularity[White, 20]], Cylinder[{{-1*10^-5, 0, 0}, {1*10^-5, 0, 0}}, r]}
 },
PlotLabel -> Chop[10^21 NIntegrate[
  Boole[cyl[{x0, y0, z0}, r, 1*10^-5] && 
    cyl[InverseFunction[xform[x, y, z, a, b]][{x0, y0, z0}], 
     1*10^-6, h]], {x0, -2*10^-5, 2*10^-5}, {y0, -2*10^-5, 
   2*10^-5}, {z0, -2*10^-5, 2*10^-5},
  Method -> {"MultidimensionalRule", "Generators" -> 9}, 
  MaxRecursion -> ControlActive[2, 4], AccuracyGoal -> 2, 
  PrecisionGoal -> 2]],
 Axes -> True,
 AxesLabel -> {"X", "Y", "Z"}],
 {{a, 0}, -1, 1},
 {{b, 1}, -1, 1},
 {{r, 1*10^-7}, 1*10^-7, 6*10^-6} ,
 {{h, 5*10^-6}, 1*10^-6, 10*10^-6},
 {{x, 0}, -1*10^-5, 1*10^-5},
 {{y, 0}, -1*10^-5, 1*10^-5},
 {{z, 0}, -1*10^-5, 1*10^-5}
]
share|improve this question
2  
What you have is a maths problem, unrelated (in its current state) to Mathematica. You need to find out how to calculate the volume of the intersection of two cylinders. Off topic here, I think. –  belisarius Nov 23 '13 at 3:20
    
Well, @MichaelE2's answer shows that my comment above was unfair –  belisarius Nov 23 '13 at 20:36

1 Answer 1

up vote 6 down vote accepted

There is a way to get Mathematica to calculate the equation of a transformed cylinder, which can then be used to calculate the volume.

First, since you're translating the cylinder, too, I rewrote your transformation to include the translation. We can also define inequalities to define the cylinder.

xform[x_, y_, z_, a_, b_] := 
 RotationTransform[a Pi/2, {1, 0, 0}] . RotationTransform[b Pi/2, {0, 1, 0}] . 
   TranslationTransform[{x, y, z}]; 
cyl[{x_, y_, z_}, r_, h_] := y^2 + z^2 < r^2 && -h < x < h;

While xform can be used to transform the Cylinder, the inverse InverseFunction[xform[x, y, z, α, β]] can be used to transform the equation of the cylinder.

In the Manipulate, you can find and display the volume of the intersection of the cylinders with

PlotLabel -> 
 Chop[10^21 NIntegrate[
    Boole[cyl[{x0, y0, z0}, r, 1*10^-5] && 
          cyl[InverseFunction[xform[x, y, z, α, β]][{x0, y0, z0}], 1*10^-6, h]],
    {x0, -h, h}, {y0, -r, r}, {z0, -r, r}, 
    Method -> {"MultidimensionalRule", "Generators" -> 9},
    MaxRecursion -> ControlActive[2, 4],
    AccuracyGoal -> 3, PrecisionGoal -> 3]]

The Chop[10^21 ...] is unnecessary, but I found it convenient. The value of the integral is quite small and NIntegrate returns a much smaller imaginary part that is only numerical error. One could use something like Chop[NIntegrate[..], 10^=30] instead. Without it, the label is too long to display nicely.

The options

MaxRecursion -> ControlActive[2, 4], AccuracyGoal -> 3, PrecisionGoal -> 3

are included because sometimes the integral is slow to converge. You can probably refine these or other NIntegrate options to suit your needs.

Update notes: When NIntegrate takes more than 5 seconds, the results are unpredictable, including messages that seem to be the result of being interrupted and not an inherent problem with the integral or the code. Alterations that improved speed:

  1. Restricted integration domain to box bounding laser beam cylinder.

  2. Added Method -> {"MultidimensionalRule", "Generators" -> 9}.

share|improve this answer
    
Thanks @Michael E2! This looks great. However, I am running into an error when I try to manipulate the rotational variables a and b in the model. NIntegrate::inumri: "The integrand Boole[y0^2+z0^2<1/100000000000000&&-(1/100000)<x0<1/100000&&1.\ (0. +Times[<<2>>]+Times[<<2>>]+Times[<<2>>])^2+1.\ (0. +Times[<<2>>]+Times[<<2>>]+Times[<<2>>])^2<1/1000000000000&&-(1/200000)<1.\ (0. +0.140901\ x0-0.31332\ y0-0.939137\ z0)<1/200000] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0,1},{0,1},{0,1}}. " Any suggestions? –  Zer0 Nov 23 '13 at 19:01
    
I've added the code I am using to implement your solution above –  Zer0 Nov 23 '13 at 19:12
    
You could try changing the domain of integration to depend on the radius r of the blue cyl. (to limit y and z) and h of the yellow (to limit x). See update. You might have to delve into the extensive NIntegrate documentation to come up with a really good solution. I run into time constraints of Manipulate (calculations should take less than 5 sec.). –  Michael E2 Nov 23 '13 at 19:46
    
@Zer0 Specified Method -- seems to help adequately. See update. –  Michael E2 Nov 23 '13 at 20:51

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.