Home

31762102746607 - ScholarWorks

1. Define whichever one of these is appropriate for the system This define needs to be the same for all the cpp files and the sxidata h file include lt fstream h gt include lt process h gt include lt stdlib h gt include lt math h gt include lt stdio h gt ifdef WIN32 include lt afxwin h gt include lt wingdi h gt endif include lt GL gl h gt include lt GL glu h gt include lt GL glaux h gt include sxi h GLOBALS GLubyte xrayData HEIGHT WIDTH GLdouble initInverse DIM DIM ifstream inFile ofstream imageOut Screen pixels NOTE this size must be the same as the dimensions of xrayData or pixel effects WILL be seen 24 int viewport 0 0 600 600 Colormaps float redmap 256 greenmap 256 bluemap 256 Viewing transformation GLdouble EYEX GLdouble 0 2 GLdouble EYEY GLdouble 0 5 GLdouble EYEZ GLdouble O 1 GLdouble CENTERX GLdouble 0 0 GLdouble CENTERY GLdouble 0 0 GLdouble CENTERZ GLdouble 0 0 GLdouble UPX GLdouble 0 0 GLdouble UPY GLdouble 0 0 GLdouble UPZ GLdouble 1 0 Scene dimensions GLfloat xMin GL 3 0 GLfloat xMax GLfloat 3 0 GLfloat yMin GLfloat 3 0 GLfloat yMax GL float 3 0 GLfloat zMin GLfloat 3 0 GLfloat 2 GL float 3 0 Mapping parameters int sensFcn 0 GLdouble exposure 1 0 char pixelType 20 MAKE IMAGE FIT
2. int MakePlasmaLoop ifstream inFile DISK testDisk0 testDisk1 testDisk2 testDisk3 char dataLine MAX DATA LINE LENGTH 1 int dataCheck int dum 0 inFile getline MAX DATA LINE LENGTH Mnitially search for the first data line It takes four lines of data to make one cylinder So read in three lines to start things off while CheckDataLine dataLine GOOD DATA linFile eof inFile getline dataLine MAX DATA LINE LENGTH sscanf dataLine amp testDisk1 radius amp testDisk1 center xVal amp testDisk1 center y Val amp testDisk1 center zVal amp testDisk1 density amp testDisk1 temperature 57 inFile getline dataLine MAX DATA LINE LENGTH while CheckDataLine dataLine GOOD DATA amp amp linFile eof inFile getline MAX DATA LINE LENGTH sscanf dataLine amp testDisk2 radius amp testDisk2 center xVal amp testDisk2 center yVal amp testDisk2 center z Val amp testDisk2 density amp testDisk2 temperature inFile getline MAX DATA LINE LENGTH while CheckDataLine dataLine GOOD DATA linFile eof inFile getline MAX DATA LINE LENGTH j sscanf dataLine amp testDisk3 radius amp testDisk3 center xVal amp testDisk3 center y Val amp testDisk3 center zVal amp testDisk3 density amp testDisk3 temperature whi
3. int main int argc char argv int statusOK SUCCESS if argc 2 AN lt lt Usage SXI lt datafile name gt lt lt endl NO FILE NAME may ultimately be used for drawing some kind of background if 1 statusOK createData xrayData else infile gt open argv 1 statusOK infile gt fail if statusOK FAILURE cout Data not available endl exit 1 statusOK SUCCESS 25 statusOK readData xrayData infile if statusOK SUCCESS myinit argv auxReshapeFunc myReshape auxMainLoop display if glGetError cout lt lt glGetError lt lt endl return statusOK Used to initialize some opengl stuff void myinit char argv ifstream inCFile sxi ini ios nocreate int i char outImage 80 GL float size Open the input data file and output image file inFile open argv 1 ios nocreate if inFile lt lt Unable to open lt lt argv 1 lt lt exiting lt lt endl exit NO DATA FILE strcpy outImage argv 1 strcat outImage bmp imageOut open outImage ios binary ios out if imageOut cerr lt lt Unable to open lt lt argv 1 lt lt bmp exiting lt lt endl exit FILE Put a light up high GLfloat pos 10 0 40 0f 100 0f 0 00 Make a substanti
4. 69 1 085 cutPI2 xParam posRootPtr gt xVal cutPI2 yParam posRootPtr gt yVal cutPI2 zParam posRootPtr gt zVal CcutPI2 offset 49 negRootPtr gt xVal cutPl1 yParam negRootPtr gt yVal cutPl1 zParam negRootPtr gt zVal cutPl1 offset pt2PI2 cutPI2 xParam negRootPtr gt xVal cutPI2 yParam negRootPtr gt cutPI2 zParam negRootPtr gt zVal cutPI2 offset If both points lie outside of the same cut plane then this ray does not intersect the cut cylinder if 21111 lt 0 22211 lt 0 241112 gt 0 pt2PI2 gt 0 goodPts FALSE j Otherwise be sure that both intersections are with portions of the cylinder that really exist after trimming by the cut planes This is accomplished by moving any point that is outside the cut plane but on the cylinder before it was cut down the ray to the point where the ray intersects the plane on the end of the cylinder else goodPts TRUE Create scrRay from the two points given as input should not matter which way we point it in or out The calculated intersection with the plane will be the same either way scrRay initPt xVal posRootPtr gt xVal scrRay initPt yVal posRootPtr gt yVal scrRay initPt zVal posRootPtr gt zVal scrRay dirVect xVal negRootPtr gt xVal posRootP
5. MONTANA STATE UNIVERSITY Synthetic x ray imager for solar plasma loops by Steven Kenneth Lundberg A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science Montana State University Copyright by Steven Kenneth Lundberg 1998 Abstract As part of developing a better understanding of the Sun models describing the magnetic phenomena responsible for the appearance of loops of plasma on the surface of the Sun are being developed The output from these models is in the form of tables of numbers Creating images from these tables that look similar to images taken of the Sun in x ray light would make using and improving the models much easier A graphical computer program has been developed to convert the numeric tables from the models into images One requirement for the program is that it be able to generate images rapidly The models predict the time development of the plasma loops and the images need to be generated and played as a movie Another requirement of the program is that it run on as many platforms as possible so that it may be used by many researchers The images generated are close approximations to the actual plasma loops described by the models They also look much like actual x ray images taken of the solar plasma loops The program is able to generate an image of several plasma loops in less than one minute This allows the generation of an entire movie consi
6. all the elements in the output matrix need to be hit so the loop will be over all elements Don t need to do anything if this element is already zero 72 if mat row DIM k 0 0 factor mat row DIM k for col 0 col lt DIM col mat col DIM factor mat col DIM row outMat col DIM factor outMat col DIM row j for row 1 lt DIM k For rows below this one Don t need to do anything if this element is already zero if mat row DIM k 0 0 factor lt mat row DIM k Matrix element on col of diag in row for col 0 col lt DIM col mat col DIM k factor mat col DIM row outMat col DIM factor outMat col DIM row End for loop over all rows Now outMat is the inverse of mat and mat is reduced to the identity matrix Copy outMat to outputMatrix for row 0 row lt DIM DIM outputMatrix row outMat row End invertMatrix Swap rows in the matrix to put the row with the largest value on the diagonal Only consider rows below the current one Performed by swapping the elements in the current matrix 3E ak ak e eak ak ae teate ae a4 ake ak ae e e ake ak ak ae ae ake ake ak ae ae ae ae 3 3e ae a4 ak ae ae ae 3e ae ak ak ae ake ae aeae ak ake ak ak ae e a ak ake e ak a e ak ak ae ae akak ak ak ak
7. 022 1 coeff z z 1 2 2 2 2 h Q33 r1 r2 zr2 2 coeffz 2 1 1 2 212 50 034 rl rl r2y zr2 const rl rl so 044 r1 r1 where zr2 15 z value where radius is r2 and at z 0 the radius is 71 All this can most easily be done in the cylinder coordinate system The screen ray is transformed to cylinder coords the intersections found and then the distance between the points is found AAA AAR HAAR REA ae ake THREE EERE AA AHH AHA AHH GLdouble CalculateIntersectionDistance int xPos int yPos GLdouble projMat DIM DIM GLdouble modelMat DIM DIM GLdouble modelMatFixed DIM DIM GLdouble botRad GLdouble topRad GLdouble height PLANE PLANE cutPI2 quadricData needs only three values since all the other 13 values of the quadric for a right truncated cylinder are either zero or one and do not change GLdouble quadricData 3 GLdouble dist 0 0 GLdouble lenVect POINTS scrPt cylStartPt cyl Vect cylEndPt 201413 posRoot negRoot 45 POINT3 tempRoot RAY int numberRoots 0 Make ray from this pixel The ray must be in cylinder coordinates but minMax and pixel position are in screen coordinates Use z near and far clipping plane values 0 0 and 1 0 Transform point on screen to point in cylinder coordinates scrPt xVal GLdouble xPos scrPt yVal GLdouble yPos scrPt zVal 0 0 Thi
8. pivot Pivot mat row if pivot row Pivot if necessary SwapRows mat pivot row SwapRows outMat pivot row Must do same thing to both matrices Make diagonal element 1 by dividing by that element Divide each element in this row in the original identity matrix outMat also if mat row DIM row 0 0 factor mat row DIM row for col 0 col lt DIM col mat col DIM row factor outMat col DIM row factor else SingularMatrix mat g EE 2 cub ML CHER DE Now zero out all elements on this column except the diagonal one by subtracting some appropriate fraction of the current row from each of the other rows Apply same operation to each of the rest of the elements in the row being operated on Also apply the same operation to each element in the original identity matrix outMat Calculate the factor to mult the ith row by to then subtract from the kth row The factor is just the value of the kth row element directly above or below the diagonal element of the ith row because the ith row element is already 1 NOTE row DIM gives the first element in a column and adding k gives the proper row below that row for k 0 k lt row k For rows above this one For all elements in the kth row The ones in the input matrix to the left of the diagonal in this row will all be zero but
9. the screen for i 0 i lt 256 i These values must be betweeen 0 and 1 SAVE THIS SET the colors look pretty good redmap i float 150 i 3 256 greenmap i float i 256 oluemap i float 0 0 290 2011 float 150 0 i 3 0 float 256 0 greenmap i float i float 256 0 bluemapf i float 0 1 redmap 0 float 0 0 greenmap 0 float 0 0 27 bluemap 0 float 0 0 glPixelMapfv GL PIXEL I 256 redmap giPixelMapfv GL PIXEL MAP I TO 256 greenmap glPixelMapfv GL PIXEL MAP I TO B 256 bluemap glPixelStorei GL UNPACK ALIGNMENT 1 glShadeModel GL FLAT glEnable GL COLOR MATERIAL glEnable GL LIGHTING glEnable GL LIGHTO glLightfv GL LIGHTO GL POSITION pos glLightfv GL LIGHTO AMBIENT amb cylQuadric gluNewQuadric Initialize the quadric for use in making cylinders and using their parameters NOTE matrix is 4 x 4 array of floats or doubles stored in column major order Thus second element in array is row 2 column 1 in matrix Can use glGet with params GL PROJECTION MATRIX amp array to get current matrix value This program intends to set up a view of a given section of the Sun say 600 000km across and 300 000km high The data in the configuration file will give the size of its view area This size will be used to scale all data values before they enter the program Currently
10. 0 0f 0 0f 1 0f blue glVertex3d xMin yMin 0 glVertex3d xMin xMax 0 glVertex3d xMax xMax 0 glVertex3d xMax yMin 0 glEnd make a small sphere at x y and z 3 glPushMatrix glTranslatef xMax 0 0f 0 0f auxWireSphere 0 02f glPopMatrix glPushMatrix glTranslatef 0 0f yMax 0 0f auxWireSphere 0 02f glPopMatrix glPushMatrix glTranslatef 0 07 0 0f auxWireSphere 0 02f glPopMatrix 55 ek IO ke ok eoe IO ICICI ok ok ok ke kc oe ok okc ke ke kc Copy elements of disk1 to disk2 void CopyDisk DISK disk1 DISK disk2 disk2 gt radius disk1 radius disk2 gt center xVal disk1 center xVal disk2 gt center yVal disk1 center yVal disk2 gt center zVal disk 1 center z Val disk2 gt density disk 1 density disk2 gt temperature disk1 temperature Check each line of data as it is read in from the input data file If data is good return the type If the data is bad exit the program with an error message int CheckDataLine char dataLine int checkedOK char revisedLine MAX
11. 1 SQRT b b 4 a c 2 a a x x b x 0 Initialize the quadric values for this cylinder quadricData 0 033 r1 r2 h 2 quadricData 1 034 043 r1 2 r1r2 h quadricData 2 044 r1 2 Where r1 botRad r2 topRad and height quadricData 0 pow botRad topRad height 2 quadricData 1 pow botRad 2 botRad topRad height quadricData 2 pow botRad 2 numberRoots GetEndPoints amp posRoot amp negRoot cylRay quadricData dist may still be zero if pts are not inside truncated cylinder The cut planes are in display coords so the points must be also if numberRoots 2 47 Transform the intersection points from cylinder coords to scene coords This is needed because they need to be checked against the clip planes and the clip plane equations are in scene coordinates Cannot use gluUnProject and gluProject here because those two functions transform cylinder coords from to screen coords not scene coordinates tempRoot TransformPoint posRoot modelMatFixed From cylinder coords posRoot tempRoot To scene coords tempRoot TransformPoint negRoot modelMatFixed From cylinder coords negRoot tempRoot scene coords if MovePointsWithinCutPlanes amp posRoot amp negRoot cutPl1 cutPI2 dist GLdouble DistBetweenPoints posRoot negRoot
12. return f static double yohkoh instr rsp 01 double temp response code for Yohkoh Al 1 int i double lgt frc static double dlog10t 0 0500002 static double tv 51 5 50000e 00 5 55000e 00 5 60000e 00 5 65000e 00 5 70000e 00 5 75000e 00 5 80000e 00 5 85000e 00 5 90000e 00 5 95000e 00 6 00000e 00 6 05000e 00 6 10000 00 6 15000 00 6 20000e 00 6 25000e 00 6 30000e 00 6 35000e 00 6 40000e 00 6 45000e 00 6 50000 00 6 55000 00 6 60000 00 6 65000 00 6 70000e 00 6 75000 00 6 80000 00 6 85000 00 6 90000 00 6 95000 00 7 00000 00 7 05000 00 7 10000 00 7 15000 00 7 20000 00 7 25000 00 7 30000 00 7 35000 00 7 40000 00 7 45000 00 7 50000 00 7 55000 00 7 60000 00 7 65000 00 7 70000 00 7 75000 00 7 80000 00 7 85000 00 7 90000 00 7 95000 00 8 00000 00 static double fv 51 2 77461 03 6 97743e 03 1 79822e 02 4 48026e 02 1 03533e 01 2 30290e 01 4 94470e 01 9 90444e 01 1 81805 00 3 11924e 00 4 96336 00 7 20120 00 9 77824 00 1 27677e 01 1 68619 01 62 2 63553e 01 5 24393e 01 9 21794e 01 132344 402 1 76815 02 2 20381e 02 2 63922 02 3 02167 02 3 28608 02 3 39777e402 3 30641 02 3 09458 02 2 83106 02 2 59439 02 2 39324e 02 2 22148 02 2 03692 02 1 87373 02 1 72428 02 1 60076 02 1 50666 02 1 43582e 02 1 38216e 02 1 34095 02 1 30807e 02
13. right xVal middle xVal lenleft lenright left x Val pl gt yParam middle yVal right yVal middle yVal lenleft lenright left yVal pl gt zParam middle zVal right zVal middle zVal lenleft lenright left zVal The plane contains the middle point pl gt offset 1 pl gt middle xVal pl gt middle yVal pl gt zParam middle zVal Return amount to extend cylinder along Z axis to allow room for cut plane to cut entirely through outside wall of cylinder The cut plane must not intersect the new end plane of the cylinder except at one point This extension will be at constant radius so that both cylinders that meet at this point will have the same cross section where they meet at the cut plane This introduces some error in the case where the cylinders are actually truncated cones The effect of the error will be to make the cylinder which has its small end at this junction a little larger than it was supposed to be and to make the cylinder which has its large end at this junction a little smaller than it was supposed to be Larger bends XS will lead to larger errors so a note will be displayed if there is more than 10 degrees bend at any junction 2777 void ExtendCyl GLfloat angZ DISK thisEnd GLfloat deltaZ angZ angZ RadToDeg Convert angZ back to
14. 00000 1 31576 2 08676 0 00000 2 08676 3 53241 4 95372 6 33112 7 64649 8 88722 10 0440 11 1068 12 0667 12 9074 13 6282 14 2237 14 6892 15 0191 15 2079 15 2623 15 1845 14 9777 14 7326 14 3923 83 2 81452 3 15818e 06 2 44063 3 39148 06 2 17009 3 59667 06 1 96982 3 77508 06 1 82274 3 92444 06 1 71365 4 04742e 06 1 63281 4 14640 06 1 57670 4 21953 06 1 54514 4 26242 06 1 52860 4 28542 06 1 52583 4 28929 06 1 53012 4 28328 06 1 54684 4 26006 06 1 58037 4 21463 06 1 63929 4 13820 06 1 72221 4 03735 06 1 83485 3 91146 06 1 98573 3 75992 06 2 19180 3 57881 06 2 47059 3 37085 06 2 85670 3 13478 06 3 42333 2 86362 06 4 28742 2 55883 06 5 70714 2 21784 06 8 36948 1 83143 06 14 5520 1 38892 06 36 5484 876406 36 5484 248285 36 5484 248285 429871 429871 1 23878 06 21 7752 1 90111 06 12 8725 2 47262 06 8 88390 2 97637 06 6 72106 3 42192 06 5 40008 3 81759 06 4 52731 4 16936 06 3 92054 4 48039 06 3 48356 4 75311 06 3 16405 4 98732 06 2 92564 5 18656 06 2 74775 5 35181 06 2 61681 5 48407 06 2 52457 5 58336 06 2 47135 5 64315 06 2 44395 5 67471 06 2 43824 5 68134 06 2 46439 5 65112 06 2 49432 5 61711 06 2 55420 5 55088 06 51 2853 51 2853 51 2853 3 57542 3 47833 3 37629 3 26971 3 16067 3 05060 2 94090 2 83292 2 72822 2 62936 2 53877 2 45946 2 39606 2 35526 2 35526 END LOOP 25 2046 26 84
15. 3 43276 3 46032 3 50339 3 61989 3 61989 END LOOP 4 80654 4 80654 4 49161 4 36181 4 26675 4 20106 4 15867 4 13266 4 11738 4 10855 4 10246 4 09577 4 08570 4 07017 4 04750 4 01626 3 97522 3 92450 3 86415 3 79451 3 73505 3 66638 0 466850 0 393197 1 28392 2 20407 3 15515 4 13185 5 12978 6 14648 7 17818 8 21885 9 05608 9 33299 10 3782 11 4185 12 4482 13 4626 14 4589 15 4323 16 3766 17 2924 18 1778 19 0261 19 8410 20 6239 21 3748 22 0990 22 8076 23 5304 22 8076 4 42275 5 48266 4 42275 3 37460 2 29320 1 15670 0 0496263 1 32993 2 68317 4 10748 5 59908 7 15612 8 76745 10 4250 12 1200 13 8437 15 5861 17 3365 19 0857 20 8250 22 1451 23 5272 1 05147 1 72146 2 41166 3 12009 3 84647 4 58671 5 33757 6 09683 6 86155 7 62815 8 24168 8 44374 9 20469 9 95867 10 7018 11 4317 12 1467 12 8441 13 5197 14 1744 14 8066 15 4109 15 9889 16 5397 17 0589 17 5435 17 9931 18 4958 17 9931 6 90290 6 74896 6 90290 6 94396 6 95393 6 92940 6 86697 6 76663 6 62968 6 45720 6 25089 6 01127 5 74349 5 45117 5 13825 4 80882 4 46750 4 11977 3 77022 3 42328 3 16214 2 89243 6 63734 7 33823 7 97835 8 55375 9 05109 9 47237 9 81835 10 0815 10 2552 10 3478 10 3627 10 3495 10 2617 10 0923 9 83369 9 49079 9 07174 8 57609 8 00040 7 35855 6 65424 5 88285 5 05716 4 18340 3 26329 2 30373 1 31576 0
16. 97 Header file for functions used to produce picture from x ray intensity input ifndef sxi sxi Must turn on one of these on the compile command line 4 WIN define UNIX gluLookAt data Also used by cylmap cpp extern GLdouble EYEX extern GLdouble EYEY extern GLdouble EYEZ extern GLdouble CENTERX extern GLdouble CENTERY extern GLdouble CENTERZ extern int viewport extern GLfloat xMin extern GLfloat extern GLfloat yMin extern GLfloat yMax extern GLfloat zMin extern GLfloat zMax extern int sensFcn extern GLdouble exposure Dimensions of synthetic xray image array in pixels NOTE WIDTH must be a multiple of four bytes or else the bmp file will not be written correctly NOTE 2 If the pixel size of the xray image buffer does not match exactly the size of the window then pixel effects WILL be seen const int WIDTH 600 const int HEIGHT 600 const int DIM 4 DIM is the dimension of the square matrix to be inverted const float PI float 3 141593 const float SQRT2 float 1 4143 const float RadToDeg float 360 0 2 0 PI 78 const int SUCCESS 1 const int FAILURE 0 const int TRUE 1 TRUE and FALSE already defined in const int FALSE 0 const int IMAGINARY INTERSECTION 101 const int LARGE ANGLE 102 const int SINGULAR MATRIX 103 const int NO DATA FILE 104 const int
17. NO BMP FILE 105 const int INVALID DATA 106 const int NO INPUT FILE NAME 107 const int MAX DATA LINE LENGTH 80 const int BLANK LINE 0 const int HEADER LINE 1 const int TITLE 2 const int END PLASMA LOOP 3 const int GOOD DATA 4 const int END OF FILE 5 extern GLUquadricObj cylQuadric Used for the construction of each cylinder extern GLubyte xrayData HEIGHT WIDTH Holds the pixel color values extern float redmap 256 Maps used to map the color values to actual colors extern float greenmap 256 extern float bluemap 256 extern GLdouble initInverse DIM DIM Needed by later function calls struct POINT3 GLdouble xVal GLdouble yVal GLdouble zVal struct DISK POINT3 center GLdouble radius GLdouble density GLdouble temperature struct PLANE By 2 0 GLdouble xParam GLdouble yParam GLdouble zParam GLdouble offset h struct RAY 79 POINT3 initPt AnitPt is starting end of RAY 20113 dirVect dirVect is the direction vector for RAY sxi cpp void myinit char argv void CALLBACK display void void CALLBACK myReshape int w int h in sxiUtils cpp int displayData GLubyte xrayData HEIGHT WIDTH void SetupAxes void void CopyDisk DISK disk1 DISK disk2 int CheckDataLine char dataLine int GoodDataLine char revisedLine int MakePlasmaLoop ifstream inFile void
18. PLANE 012 double numer denom GLfloat dot numer pl1 xParam pl2 xParam pll yParam pl2 yParam pll zParam pl2 zParam denom sqrt pow pl1 xParam 2 pow pll yParam 2 pow pl1 zParam 2 sqrt pow pl2 xParam 2 pow pl2 yParam 2 pow pl2 zParam 2 dot GLfloat numer denom return dot Make planar equation from PLANE set up direction void SetEqn PLANE GLdouble eqn GLint up eqn 0 GLdouble up pl xParam eqn 1 GLdouble up pl yParam eqn 2 GLdouble up pl zParam eqn 3 GLdouble up pl offset Calculate the distance between two points GLfloat DistBetweenPoints 201113 pt1 POINT3 pt2 38 return GLfloat sqrt pow pt2 xVal pt1 xVal 2 pow pt2 yVal ptl yVal 2 pow pt2 zVal pt1 zVal 2 39 APPENDIX C CYLMAP CPP 40 cylmap cpp Steve Lundberg 7 20 97 This file contains the functions used to perform ray tracing from the screen image of a cylinder in the position it is initially created z axis symetric bottom on origin through the cylinder The length of the ray through the cylinder is then used to scale the color of the pixel where the ray originated include lt stdlib h gt include lt fstream h gt include lt process h gt include lt math h gt ifdef WIN32
19. a ak ake ak ake ak ake ake 29 Called when the window is first opened and whenever the window is reconfigured void CALLBACK display int statusOK i j GLdouble initMat DIM DIM int dataCheck ffstream inFile xdata in ios nocreate glClear GL_COLOR_BUFFER_BIT Takes care of background painting Initialize xrayData to black for i 0 lt HEIGHT for j 0 j lt WIDTH j xrayData i j GLubyte 0 Place a color table in bottom left of pixel array xrayData for i 0 i lt 20 for j 0 j lt 256 xrayData i j GLubyte j SetupAxes glPolygonMode GL FRONT AND BACK GL FILL Get the current MODELVIEW matrix and store it s inverse for later use glGetDoublev MODELVIEW MATRIX InvertMatrix initMat initInverse 52 05 ON CES ME PEDI PES VE NE This section loops as long as there are still records in the input file The first three records enough to set up the following loop were read in the init function Now we repeatedly move each disk over one and read in one new record Using the four records disks we construct one cylinder At the end of the file we are done Here we read in a line from the input data file as a string Then check that line to see if it is valid data or if it is the END OF LOOP code Or some other code If the line i
20. ak ak ak aeae akak ak void SwapRows GLdouble int pivot int row int col GLdouble temp for col 0 col lt DIM col temp mat pivot col DIM mat pivot col DIM mat row col DIM mat row col DIM temp 73 Exit if the matrix is singular void SingularMatrix GLdouble mat cout lt lt endl lt lt Matrix is singular Cannot continue lt lt endl PrintMatrix mat exit SINGULAR MATRIX Hehe he ok Find the largest element below or on the diagonal indexed by row The pivot must be chosen at the time it is needed from those available int Pivot GLdouble mat int row int pivot ckRow GLdouble max if row lt DIM 1 The last row cannot swapped max fabs mat row DIM row Set diagonal element as the max pivot row for ckRow rowt 1 ckRow lt DIM ckRow Check rows below this row if fabs mat row DIM ckRow gt max max fabs mat row DIM ckRow pivot ckRow else H pivot DIM 1 The last row cannot swapped with anyone return pivot Print the given matrix void PrintMatrix GLdouble mat int row col for row 0 row lt DIM row for col 0 col lt DIM col printf 4 2f mat row col D
21. bfSize 4 memcpy buff 6 amp bitMapFileHeader bfReserved1 2 memcpy buff 8 amp bitMapFileHeader bfReserved2 2 memcpy buff 10 amp bitMapFileHeader bfOffBytes 4 imageOut write buff 14 Put the bmiHeader in the bmp file memopy buff amp bmiHeader biSize 4 memcpy buff 4 amp bmiHeader biWidth 4 memcpy buff 8 amp bmiHeader biHeight 4 59 12 amp bmiHeader biPlanes 2 14 amp bmiHeader biBitCount 2 memcpy buff 16 amp bmiHeader biCompression 4 memcpy buff 20 amp bmiHeader biSizeImage 4 memcpy buff 24 amp bmiHeader biXPelsPerMeter 4 memcpy buff 24 amp bmiHeader biYPelsPerMeter 4 28 amp bmiHeader biClrUsed 4 memcpy buff 32 amp bmiHeader biClrImportant 4 imageOut write buff 40 Put the colormap in the bmp file color 3 0x00 The fourth byte of the quad is zero for int i 0 i lt 256 i color 0 unsigned char bluemap i 255 0 color 1 unsigned char greenmap i 255 0 color 2 unsigned char redmap i 255 0 imageOut write color 4 Put the pixels in the bmp file Starting with the lower left pixel and working one row at a time The length of a row must be a multiple of 4 bytes for i 0 i lt HEIGHT i for int j 0 j lt WIDTH j imageOut write amp xrayData i j 1 j imageOut close Map the temp
22. has an associated radius a density factor and a temperature factor By connecting the points with an envelope that extends a radius from each point a three dimensional loop is generated Since the plasma density and temperature typically vary au the length of the plasma loop these properties are also associated with each point SXI averages the density and temperature attributes at each two connected points to obtain the values used for that portion of the loop The camera in the satellite uses different exposure times and filters to effectively image different parts of the solar plasma loops SXI uses the measured properties of the filters in Yohkoh to create images that embody the same responses to x ray energy Yohkoh images the Sun with a CCD camera This yields pictures with a fixed 5 pixel resolution SXI is designed so that the pixel size in the images can be set to any desired size or allowed to adjust so that the entire plasma loop is within the image One primary aspect of the images taken by Yohkoh and modeled by the solar magnetic models is the time development of the plasma loops Through interactions with other loops and under control of the transport phenomena inside the Sun the plasma loops grow merge move and change temperatures This dynamic aspect of the plasma loops is particularly important for verifying the model accuracy Thus SXI is designed to be able to generate images quickly enough to create a movie wi
23. if dist 0 8 DEBUG CLIP PLANE cout large dist dist endl j return dist ARR ke RH e joke kk je ke ke oe ke ke ke e ok ke ke ke ke ke ke e ok ke ke e je ke ke je je ke ke oe oko ke KAA RAK RAHA HEHE KEK First calculate the K values then set up the solutions using the quadratic equation b 8 2 4ac 2a or K1 SQRT K1 2 4 K2 K0 2 K2 int GetEndPoints POINT3 posRootPtr POINT3 negRootPtr RAY cylRay GLdouble quadricData GLdouble K0 K1 K2 detmnt posRootT negRootT int numRoots 0 pow cylRay initPt xVal 2 pow cylRay initPt yVal 2 pow cylRay initPt zVal 2 quadricData 0 quadricData 2 2 cylRay initPt zVal quadricData 1 K1 2 cylRay dirVect xVal cylRay initPt xVal 2 cylRay dirVect yVal cylRay initPt yVal 2 cylRay dirVect zVal cylRay initPt zVal quadricData 0 quadricData 1 K2 pow cylRay dirVect xVal 2 pow cylRay dirVect yVal 2 pow cylRay dirVect zVal 2 quadricData 0 detmnt pow K1 2 4 K2 if detmnt 0 amp amp 0 if detmnt 0 48 numRoots 1 tangential intersection distance 0 else numRoots 2 posRootT K1 sqrt detmnt 2 K2 negRootT K1 sqrt detmnt 2 K2 posRootPtr gt xVal cylRay initPt xVal GLfloat posRootT cylRay dirVec
24. include lt afxwin h gt include lt wingdi h gt endif include lt GL gl h gt include lt GL glu h gt include lt GL glaux h gt include sxi h DEBUG GLOBAL static int DEBUG MOVE POINTS 0 static int DEBUG CLIP PLANE 0 First calc direction vector for screen into scene Since we are using parallel projection there is only one direction vector for the entire scene Second calc bounding rectangle for cylinder in screen coords Third for each pixel in bounding rectangle determine if it intersects with the cylinder This requires mapping the screen pixels and rays to the cylinder coordinates then calculating the intersections Fourth if the cylinder intersection is outside the clip planes for the ends of this cylinder replace this intersection with the intersection of the ray with the clip plane But if the second ray intersection with the cylinder is also outside the same clip plane then this ray doesn t intersect the clipped cylinder at all Fifth calculate the distance between the two points of intersection Sixth set the xray image buffer pixel color proportional to the distance This function is called with the modelMatrix for the construction of the cylinder i e cylinder axis is z axis bottom of cylinder at z 0 k
25. int wpos 1 lt 2 j if xrayData i j gt GLubyte 220 count return count void increaseLocal9 GLubyte xrayData HEIGHT WIDTH int hpos int wpos int amt m temp for int hpos 1 i lt 2 i for int j wpos 1 lt 2 j int xrayData i j temp amt if temp gt 255 temp 255 else if temp 0 temp 0 67 xrayData i j GLubyte temp int readData GLubyte xrayData HEIGHT WIDTH ifstream infile for int i 0 i lt HEIGHT i for int j 0 j lt WIDTH j return infile gt fail infile gt get xrayData j i Pug 68 APPENDIX E MATRIX CPP fe matrix cpp Steve Lundberg 11 12 97 This file contains utility functions having to do with matrix operations Functions included InvertMatrix uses Gauss Jordan elimination method SwapRows used by InvertMatrix SingularMatrix used by InvertMatrix Pivot used by InvertMatrix PrintMatrix MatrixMult TransformPoint gt include lt stdlib h gt include lt fstream h gt include lt process h gt include lt math h gt ifdef WIN32 include lt afxwin h gt include lt wingdi h gt endif include lt GL gl h gt include lt GL glu h gt include lt GL glaux h gt include sxi h Invert matrix calc the inverse of the current model
26. radians 8 must be lt 45 deg deltaZ GLfloat thisEnd radius GLfloat tan angZ Check the angle between two cylinders If it is greater than 4506 exit the program with an error If it is larger than 10068 display a warning void CheckAngle GLfloat angZ POINT3 pt if angZ gt 45 0 For flat end the cut plane angle 15 0 deg cerr ERROR the angle between cylinders is too great endl 37 lt lt The angle is lt lt angZ lt lt degrees lt lt endl cerr lt lt The center point is at x y z lt lt pt xVal lt lt lt lt pt yVal lt lt lt lt pt zVal lt lt endl exit LARGE_ANGLE j if angZ 10 0 Leads to larger errors notify user cout lt lt NOTE the angle between adjacent cylinders lt lt endl cout lt lt is greater than 10 degrees This may yield lt lt endl cout lt lt unacceptably large errors in the mapping process lt lt endl cout lt lt The center point is at x y z lt lt pt xVal lt lt lt lt pt yVal lt lt lt lt pt zVal lt lt endl cout lt lt The angle is lt lt angZ lt lt endl 9C He ke oko ke ke ok ok ke ke koe ke kc kc ke e ke ke ke ke ke ke ok oe ke ok oe ke je je Calculate the dot product for two vectors defined by two planes k okk kkk kkk GLfloat DotProd PLANE
27. this scaling is not implemented The internal representation in the program will be of a volume that is 6 units across X and Y and 3 units high Z Data points that are outside of this area will be flagged and printed out This range checking will be done as the data are read the first time HE AE eH AE He ee ee ee ade ste ee ee ee ee ee ee A ade ake ae sle 7 Called when the window is first opened and whenever the window is reconfigured void CALLBACK myReshape int w int giMatrixMode PROJECTION glLoadIdentity initialize the matrix GLfloat xOrthoMin xOrthoMax yOrthoMin yOrthoMax zOrthoMin zOrthoMax GL float adj Scale the viewing volume so the entire x y plane will fit When viewing at an angle to the x axis or y axis that is not 90 or 0 degrees then viewing volume needs to be expanded so that the corners of the x y plane do not get cut off adj GLfloat 1 0 EYEX GLfloat 0 0 EYEY GLfloat 0 0 28 if fabs EYEX gt fabs EYEY adj GLfloat 1 0 0 4 fabs EYEY EYEX else adj GLfloat 1 0 0 4 fabs EYEX EYEY if EYEZ GLfloat 0 0 adj adj GLfloat fabs EYEZ zMax xOrthoMin adj xMin xOrthoMax adj xMax yOrthoMin adj yMin yOrthoMax adj yMax zOrthoMin adj zMin zOrthoMax adj zMax define the viewing volume glO
28. 0 6 05000 00 6 10000 00 6 15000e 00 6 20000e 00 61 6 25000e 00 6 30000 00 6 35000 00 6 40000 00 6 45000e 00 6 50000e 00 6 55000 00 6 60000e 00 6 65000 00 6 70000 00 6 75000e 00 6 80000 00 6 85000e 00 6 90000e 00 6 95000e 00 7 00000e 00 7 05000e 00 7 10000 00 7 15000e 00 7 20000e 00 7 25000e 00 7 30000e 00 7 35000e 00 7 40000e 00 7 45000e 00 7 50000e 00 7 55000e 00 7 60000e 00 7 65000e 00 7 70000e 00 7 75000e 00 7 80000e 00 7 85000e 00 7 90000e 00 7 95000e 00 8 00000e 00 static double fv 51 1 00895e 02 2 40323e 02 5 79735e 02 1 35350e 01 2 93222e 01 6 12964e 01 1 24098e 00 2 35314e 00 4 12038e 00 6 77274e 00 1 03732 01 1 44994e 01 1 89599 01 2 37722 01 2 98577 01 4 26198 01 7 66350 01 1 27973 02 1 79364 02 2 35745 02 2 90719 02 3 45278 02 3 92810 02 4 25096 02 4 37849 02 4 24666 02 3 96204 02 3 61350 02 3 30211 02 3 03990 02 2 81867 02 2 58477 02 2 37962 02 2 19285 02 2 03900 02 1 92195 02 1 83391 02 1 76720 02 1 71596 02 1 67501 02 1 64061 02 1 60991 02 1 58136 02 1 55410 02 1 52718 02 1 50017 02 1 47336 02 1 44690 02 1 42077 02 1 39497 02 1 369636 02 Igt log10 temp i int floor Igt tv 0 dlog10t if 1 lt 0 return 0 0 if 1 gt 49 return fv 50 frc Igt tv i tv i 1 tv i f frc fv i 1 1 0 frc fv i
29. 00 6 25000e 00 6 30000e 00 6 35000e 00 6 40000e 00 6 45000e 00 6 50000e 00 6 55000e 00 6 60000e 00 6 65000e 00 6 70000e 00 6 75000e 00 6 80000e 00 6 85000e 00 6 90000e 00 6 95000e 00 7 00000e 00 7 05000e 00 7 10000e 00 7 15000e 00 7 20000e 00 7 25000e 00 7 30000e 00 7 35000e 00 7 40000e 00 7 45000e 00 7 50000e 00 7 55000e 00 7 60000e 00 7 65000e 00 7 70000e 00 7 75000e 00 7 80000e 00 7 85000e 00 7 90000e 00 7 95000e 00 8 00000e 00 static double fv 51 1 28114e 12 1 38595e 11 4 80228e 10 4 41040e 09 4 73539e 08 3 58486e 07 2 96677e 06 1 56400e 05 7 17745e 05 2 93496e 04 1 09007e 03 3 56551e 03 1 01255e 02 2 59010e 02 6 25827e 02 1 62521e 01 4 43389e 01 9 62890e 01 1 68332e 00 2 67787 00 3 61954 00 4 81358e 00 6 08095e 00 7 33778 00 8 59615e 00 9 85502 00 1 12021 01 1 24850e 01 1 35449 01 1 41905 01 1 45598 01 1 45653 01 1 43720 01 1 38685 01 1 32029 01 1 25454 01 1 19659 01 1 14821 01 1 10789 01 1 07338 01 1 04341 01 1 01662 01 9 91788 00 9 68007 00 9 45109 00 9 22976 00 9 01683 00 8 81230 00 8 61912 00 8 43765 00 8 26690 00 Igt log10 temp i int floor Igt tv 0 dlog10t if i lt 0 return 0 0 1 i gt 49 return fv 50 frc Igt tv i tv i 1 f fre fv i 1 1 0 frc fv i return f static double yohkoh_instr_rsp_05 double temp response code
30. 1 28051e 02 1 25597e 02 1 23322e 02 1 21157e 02 1 19026e 02 1 16894 02 1 14783 02 1 12705e 02 1 10659 02 1 08644e 02 1 06669 02 Igt log10 temp i int floor Igt tv 0 dlog10t 1 i lt 0 return 0 0 if 1 gt 49 return fv 50 frc Igt tv i tv i 1 tv i f fre fv i 1 1 0 frc fv i return f static double yohkoh instr rsp 02 double temp response code for Yohkoh Al Mg int 1 double Igt frc static double dlogl0t 0 0500002 static double tv 51 5 50000 00 5 55000e 00 5 60000 00 5 65000 00 5 70000e 00 5 75000e 00 5 80000e 00 5 85000 00 5 90000 00 5 95000e 00 6 00000e 00 6 05000e 00 6 10000e 00 6 15000e 00 6 20000e 00 6 25000e 00 6 30000e 00 6 35000e 00 6 40000 00 6 45000e 00 6 50000 00 6 55000 00 6 60000 00 6 65000 00 6 70000e 00 6 75000 00 6 80000 00 6 85000 00 6 90000 00 6 95000 00 7 00000 00 7 05000 00 7 10000 00 7 15000 00 7 20000 00 7 25000 00 7 30000 00 7 35000 00 7 40000 00 7 45000 00 7 50000 00 7 55000 00 7 60000 00 7 65000 00 7 70000 00 7 75000 00 7 80000 00 7 85000 00 7 90000 00 7 95000 00 8 00000e 00 static double fv 51 2 07817e 04 6 80402e 04 2 05662e 03 5 86829e 03 1 55755e 02 3 94193e 02 9 45655e 02 2 09322e 01 4 18741e 01 7 75536e 01 1 32438 00 2 06871 00 3 04000e 00 4 31498 00 6 2
31. 1259 00 1 08623 01 2 38189 01 4 39982 01 6 52788 01 8 98284 01 1 14151 02 1 39553 02 1 62925 02 1 80595 02 1 90506 02 1 89800 02 1 82796 02 1 72811 02 1 63677 02 1 55329 02 1 47332e 02 1 37159 02 1 27405 02 1 17853 02 1 096466 02 1 03261 02 9 84019 01 9 47082 01 9 18763 01 8 96294 01 8 77615 01 8 61131 01 8 45952 01 8 31563 01 8 17426 01 8 03298 01 7 89331 01 7 75595 01 7 62076 01 7 48760 01 7 357096 01 lgt log10 temp i int floor igt tv 0 107 if i lt 0 return 0 0 63 if 1 gt 49 return fv 50 frc Igt tv i tv i 1 tv i f fre fv i 1 1 0 frc fv i return f static double yohkoh instr rsp 03 double temp response code for Yohkoh Mg3 int i double Igt frc static double dlog10t 0 0500002 static double tv 51 5 50000e 00 5 55000 00 5 60000e 00 5 65000 00 5 70000e 00 5 75000e 00 5 80000e 00 5 85000e 00 5 90000e 00 5 95000e 00 6 00000e 00 6 05000e 00 6 10000e 00 6 15000e 00 6 20000e 00 6 25000e 00 6 30000 00 6 35000 00 6 40000e 00 6 45000e 00 6 50000e 00 6 55000e 00 6 60000 00 6 65000e 00 6 70000e 00 6 75000e 00 6 80000 00 6 85000 00 6 90000e 00 6 95000e 00 7 00000e 00 7 05000 00 7 10000 00 7 15000 00 7 20000e 00 7 25000 00 7 30000 00 7 35000 00 7 40000 00 7 45000 00 7
32. 2 z0 z0 Q33 044 2 z0 Q34 Where The Q s are the non zero elements of a 4x4 matrix the quadric form of the circular truncated cone 46 ax ay and az are the direction vector elements for the ray X0 0 and 20 are the starting end point values for the ray ref Graphics Gems III p 280 The general form for a quadric surface is F x y z ax 2 2bxy 2cxz 2dx 2 2fxy 2gy hz 2 22 0 The matrix form of a quadric is ja be befg y 1 2 y 2 1 cfhi z XQX 0 a 012 021 b 014 041 d etc For a right truncated cone with z the axis of symmetry the values of the non zero elements Q11 022 1 033 r1 r2 2 h 2 034 043 r1 2 r1r2 h and 044 r1 2 Where rl is the radius at z 0 12 is the radius 2 h and is height of the cylinder NOTE While it should be possible to derive the quadric for the cylinder with ends cut by arbitrary planes since I have not done that I will simply check all points to be sure they are not outside of the end planes That means that some intersections will be with the end planes rather than with the cylinder itself Solve for the values of t using the quadratic equation Must first make sure that the portion in the square root is positive Must also make sure that the intersections lie both within the cylinder volume and within the cut plane bounds of the cylinder
33. 50000 00 7 55000 00 7 60000 00 7 65000 00 7 70000 00 7 75000 00 7 80000 00 7 85000 00 7 90000 00 7 95000 00 8 00000 00 static double fv 51 5 41356e 06 2 20848e 05 8 12592e 05 2 78793e 04 8 90374e 04 2 62938e 03 7 21077e 03 1 81674e 02 4 22954e 02 9 29136e 02 1 93202e 01 3 78956e 01 7 05772e 01 1 27351e 00 2 36600e 00 5 83541e 00 1 65382e 01 3 36555e 01 5 19335e 01 7 325686 01 9 46421e 01 1 16887e 02 1 37194e 02 1 52145e 02 1 59679e 02 1 57264e 02 1 48831e 02 1 37717e 02 1 27556e 02 1 18359e 02 1 09693 02 9 94492e 01 8 99840 01 8 11926 01 7 39861 01 6 86123 01 6 46433 01 6 16809 01 5 94385 01 5 76801 01 5 62177 01 5 49208 01 5 37420 01 5 26525 01 5 16035 01 5 05703 01 4 95682 01 4 86035 01 4 76696 01 4 67631 01 4 58902 01 lgt log10 temp i int floor Igt tv 0 dlog10t if i lt 0 return 0 0 if i gt 49 return fv 50 frc Igt tv i tv i 1 f 1 0 frc fv i return f static double yohkoh instr rsp 04 double temp response code for Yohkoh 112 int i double lgt frc static double dlog10t 0 0500002 64 static double tv 51 5 50000 00 5 55000e 00 5 60000e 00 5 65000e 00 5 70000e 00 5 75000e 00 5 80000e 00 5 85000e 00 5 90000e 00 5 95000e 00 6 00000e 00 6 05000e 00 6 10000e 00 6 15000e 00 6 20000e
34. 81 28 4537 30 0155 31 5340 33 0083 34 4384 35 8238 37 1646 38 4644 39 7252 40 9492 42 1372 43 2919 42 1372 2 57179 2 26702 1 98051 1 71547 1 47207 1 25081 1 05186 0 875258 0 720907 0 586496 0 470527 0 371190 0 286779 0 214157 0 286779 13 8711 13 2430 12 5154 11 6909 10 7830 9 79986 8 74872 7 63563 6 46584 5 24841 3 98873 2 69172 1 36063 5 12227e 07 1 36063 84 2 65903 5 44036 06 2 80522 5 29671 06 3 00174 5 12038106 3 26786 4 90747 06 3 62582 4 658936 06 4 11458 4 37347e 06 4 80019 4 0491 1 06 5 80406 3 68234e 06 7 37048 3 26769 06 10 0246 2 801926 06 15 1810 2 27688e 06 27 9249 1 67878e 06 83 1051 973140 83 1051 146961 83 1051 146961
35. 9 01 1 03053 01 1 03637 01 1 03947 01 1 039936 01 Igt log10 temp i int floor Igt tv 0 dlog10t if i lt 0 return 0 0 if 1 gt 49 return fv 50 frc Igt tv i tv i 1 tv i f fre fv i 1 1 0 frc fv i return f The functions below not used but may be of some value in later attempts manipulate the pixel buffers int createData GLubyte xrayData HEIGHT WIDTH srand 97 Fixed seed value results in same set every time for int i 0 i lt HEIGHT i for int j 0 j lt WIDTH j lye xrayData i j GLubyte rand 96255 xrayData i j GLubyte 0 xrayData i j GLubyte j This shows the total range of colors commented out below for j 100 j lt 150 j xrayData i j GLubyte 150 fot 150 j lt 200 j xrayData i j GLubyte 200 for 200 j lt WIDTH j xrayData i j GLubyte 255 end of commented out 66 return 1 void playData GLubyte xrayData HEIGHT WIDTH for int i 1 1 lt 99 i 3 for int j 1 j lt WIDTH 1 j 3 if num checkLocal9 xrayData i gt 4 increaseLocal9 xrayData i j 50 else if num lt 3 increaseLocal9 xrayData i j 50 int checkLocal9 GLubyte xrayData HEIGHT WIDTH int hpos int wpos int count 0 for int i hpos 1 i lt 2 i for
36. DATA LINE LENGTH 1 J Strip any leading white space sscanf t n s amp revisedLine sscanf dataLine 0 amp revisedLine lt lt revLine lt lt revisedLine lt lt lt lt endl strcpy revisedLine dataLine the line had nothing but whitespace get the next line if strlen revisedLine 0 checkedOK BLANK LINE else if revisedLine 0 H checkedOK HEADER LINE else revisedLine 0 T checkedOK TITLE else if revisedLine 0 checkedOK END PLASMA LOOP else if GoodDataLine revisedLine checkedOK GOOD DATA else 56 lt lt nIncorrect line in input lt lt endl cerr dataLine endl exit INVALID DATA return checkedOK Determine that all six data values are within the required ranges And that there are six values TAIRA Ekk int GoodDataLine char revisedLine int checkedOK FALSE num double rad x y z den temp num sscanf revisedLine amp rad amp x amp z amp den amp temp if 0 0 lt rad amp amp 0 0 den amp amp 0 0 lt temp amp amp temp lt 15 6 amp amp num 6 checkedOK TRUE return checkedOK Make a plasma loop from individual tapered truncated cylinders
37. IM 74 j cout endl j j Multiply inMat1 inMat2 gt outMat A TT TA AAA AAA AA AAA AAT void MatrixMult GLdouble inMat1 GLdouble inMat2 GLdouble outMat int row col k for row 0 row lt DIM row for col 0 col lt DIM col outMat row col DIM 0 for k 0 k lt DIM k outMat row col DIM inMat1 row k DIM inMat2 col DIM k Function to mult matrix by POINT3 Set w 1 for the point Definitions of openGL require this to be done as M Pt not Pt M POINT3 TransformPoint POINT3 pt GLdouble matrix GLdouble inPt DIM outPt DIM int row col POINT 3 tranPt Initialize inPt from pt inPt 0 pt xVal inPt 1 pt yVal inPt 2 pt zVal inPt 3 1 0 for row 0 row lt DIM row outPt row 0 0 for col 0 col lt DIM col outPt row matrix row col DIM inPt col Copy to tranPt and return it tranPt xVal GL float outPt 0 tranPt yVal GLfloat outPt 1 tranPt zVal GLfloat outPt 2 if outPt 3 1 0 cout lt lt w not 1 0 for this Xform lt lt endl PrintMatrix matrix return tranPt 76 APPENDIX F SXLH 77 sxi h Steve Lundberg 2 15
38. R CPP 31 APPENDIX C als hg acd es kee eia oie is 39 APPENDIX D APPENDIX E MATRIX CPP 68 APPENDIX F AR P dese Boca E eiat 76 APPENDIX G EXAMPLE INITIALIZATION AND DATA FILES 81 vi LIST OF FIGURES Figure 1 Tapered Cylinder Building Block ead Figure 2 Creating Bounding Box for Cylinder 11 Figure 3 Ray Intersection with Cut Plane 12 Figure 4 Grayscale Version of Synthetic X Ray 14 vil ABSTRACT As part of developing a better understanding of the Sun models describing the magnetic phenomena responsible for the appearance of loops of plasma on the surface of the Sun are being developed The output from these models is in the form of tables of numbers Creating images from these tables that look similar to images taken of the Sun in x ray light would make using and improving the models much easier graphical computer program has been developed to convert the numeric tables from the models into images One requirement for the program is that it be able to generate images rapidly The models predict the time development of the plasma loops and the images need to be generated and played as a movie A
39. SaveImage GLubyte xrayDataHEIGHT WIDTH ofstream imageOut void ApplySensitivityFcn int fcn GLdouble temp in cylinder cpp void MakeCylinder DISK prior DISK bot DISK top DISK after GLfloat DistBetweenPoints POINT3 pt1 POINT3 pt2 void MoveToOrigin POINT3 pt void AlignZAxis POINT3 botPt POINT3 topPt void CutPlane POINT3 left POINT3 middle POINT3 right PLANE pl void ExtendCyl GLfloat angZ DISK thisEnd GLfloat deltaZ void CheckAngle GLfloat angZ POINT3 GLfloat DotProd PLANE pl1 PLANE pl2 void SetEqn PLANE pl GLdouble eqn GLint up GLfloat DistBetweenPoints POINT3 pt1 POINT3 pt2 matrix cpp void InvertMatrix GLdouble inputMatrix 16 GLdouble outputMatrix void SwapRows GLdouble mat int pivot int i void SingularMatrix GLdouble mat int Pivot GLdouble mat int row void PrintMatrix GLdouble mat void MatrixMult GLdouble inMat1 GLdouble inMat2 GLdouble outMat POINTS TransformPoint POINT3 cyIPt GLdouble Matrix in cylmap cpp void MapCylinderPixels GLdouble botRad GLdouble topRad GLdouble height PLANE PLANE cutPI2 GLdouble density GLdouble sensitivity1 GLdouble sensitivity2 void MinMax GLdouble minMax POINT3 scrPt int init GLdouble CalculateIntersectionDistance int xPos int yPos GLdouble projMat GLdouble modelMat GLdouble modelMatFixed GLdouble botRad GLdouble topR
40. The actual picture of the loop drawn on the screen is used only for pixel region determination it is not used for coloring the synthetic xray image that is all done with the mathematical descriptions of the cylinders 2 The openGL projection matrix used to map the current cylinder to the screen can be inverted to map the ray from screen coordinates to cylinder coordinates where the mathematics for the cylinder are simpler 3 Determining which screen pixels contain information from the cylinder 23 can be done by drawing a rectangle around the cylinder in screen coordinates All pixels inside this rectangle will be used to start rays into the scene Each of these rays will be transformed using the inverse of the projection matrix and then intersected with the cylinder 4 If the ray does intersect the cylinder The intersection must also be between the two cylinder end planes the distance between the two intersections will be calculated The distance is further refined by a density factor and a temperature factor So distances through regions where the density and or temperature are low will not be as bright as regions where the temperature and density are relatively high 5 The depth of penetration information will be converted to an appropriate value and added to the current pixel value Thus if more than one cylinder is penetrated by the same ray the sum of the distances through all cylinders will determine the pixel color
41. The temperature variable is x ray luminosity convolved with the response function of the filter and recorder being simulated X rays from plasma of different temperatures have different characteristic wavelengths and each filter material has a different response to x ray wavelength The temperature variable is calculated by using the temperature in degrees Kelvin read from the input file and a function built in to the program The function calculates the expected number of counts in an individual detector bin from plasma of the given temperature and normalized for exposure time electron number da depth of the plasma path and the area on the Sun that is imaged by the specific detector bin When this value is multiplied by the actual exposure time electron number density and plasma path depth it gives the expected counts per pixel for the synthetic x ray image One further adjustment needs to be made The Yohkoh camera uses a charge coupled device CCD to collect the x ray photons and the counts are digitized to 12 bit accuracy This program is using 8 bit values to store the pixel colors To make the SXI images saturate at the same x ray intensity as the Yohkoh camera the counts per pixel are divided by 16 to allow for the greater dynamic range of Yohkoh Mapping Transformations The mapping process requires both coordinate transformations and intersection calculations The coordinate systems will be discussed first The final synthe
42. ad GLdouble height PLANE PLANE cutPI2 int GetEndPoints POINT3 posRootPtr POINT3 negRootPtr RAY cylRay 80 GLdouble quadricData int MovePointsWithinCutPlanes POINT3 posRootPtr POINT3 negRootPtr PLANE cutPl1 PLANE cutPI2 void IntersectRayWithCutPlane POINT3 rootPtr PLANE RAY cylRay void AdjustPixelColor int xPos int yPos int dist GLdouble GLdouble temp int tempFcn endif 81 APPENDIX G EXAMPLE INITIALIZATION AND DATA FILES 82 File sxi ini 0 0 600 600 0 0 0 200 20000001 1 08 FIXED PIXEL SIZE 0 14 1500 File format line 1 windowX windowY windowWidth windowHeight line 2 EYEX EYEY EYEZ CENTERX CENTERY CENTERZ UPX UPY UPZ line 3 Sensitivity function index exposure value pixel sizing size line 4 xOffset yOffset zOffset Line 1 controls the size in pixels of the window on the screen This needs to be the same size as the final synthetic xray image will be or pixel artifact effects will be seen Line 2 controls the viewpoint from which the scene is drawn The three EYE coordinates and the three CENTER coordinates determine the viewing angle The UP coordinates determine which direction is shown as up on the screen Line 3 sets the sensitivity function to be used This relates to the type of filter in use And sets the exposure which relates to the brightness ofthe image And defines the type of scalin
43. added in The temperature is scaled proportional to the distance along the path between the cylinder end points An exposure time factor is also used to allow for enhancing different sections of the image This section is heavily openGL dependent whereas most of the preceeding work is not The distance values are scaled to a color value and written to a pre defined storage buffer xPos and yPos are pixel numbers First make sure the bounding rectangle is entirely within the displayable image Trim it to fit if it isn t if minMax 0 0 43 minMax 0 0 vt minMax 2 0 0 if minMax 1 gt WIDTH ON WIDTH rl minMax 3 gt HEIGHT minMax 3 HEIGHT for yPos int minMax 2 yPos lt int minMax 3 yPos for xPos int minMax 0 xPos lt 1 dist CalculateIntersectionDistance xPos 4105 projMat modelMat modelMatFixed botRad topRad height 69 1 cutPI2 Scale distance by density and sensor sensitivity The distance is 1 0 for an intersection distance 1 6 of horizontal span Density is the square of the electron number density read from the input file The exposure variable allows for brightness adjustment The result is divided by 16 to make the count match what would be stored in a 12 bit value This value is 8 bits counts int dist xMax xMin 6 density sensitivi
44. agnetic models It produces synthetic xray pictures from the size density and temperature of the loops This program is capable of rendering several plasma loops in one image Although this code is written to be compiled with a compiler nearly all of the code is standard C olf Process A Use the points defined in the input data together with the radii also defined in the input data to produce tapered cylinders with ends cut at an angle that neatly meets the next cylinder B From the view point of the screen view for each screen pixel that is in a cylinder Color that pixel proportional to the distance through the cylinder in the direction from the screen pixel perpendicular to the screen Another way to say this is that the pixel is colored proportional to the distance a ray from the screen pixel and perpendicular to the screen will travel from the point where it enters the cylinder to the point where it leaves the cylinder Since all cylinders are convex objects no ray can go through more than one place on an individual cylinder It is possible however for a given ray to penetrate more than one cylinder This simply means that the pixel color will be determined by the total depth of all cylinders that ray passes through 1 These rays are mathematically determined by the coordinates of the screen pixel and the look at direction Since this is an orthogonal projection all rays go in the same direction
45. al glow from all directions GLfloat amb 0 6f 0 6f 0 6f 1 00 Read in the configuration values if inCFile cerr lt lt Unable to open sxicfg in continuing by using lt lt endl lt lt default values for configuration variables lt lt endl 26 else inCFile gt gt viewport 0 gt gt viewport 1 gt gt viewport 2 gt gt viewport 3 inCFile gt gt EYEX gt gt EYEY gt gt EYEZ gt gt CENTERX gt gt CENTERY gt gt CENTERZ gt gt UPX gt gt UPY gt gt UPZ inCFile gt gt sensFcn gt gt exposure gt gt pixelType gt gt size pixelType 19 0 if pixelType FIXED AREA SIZE 15 FIXED PIXEL SIZE so calculate area xMax lt size WIDTH float 2 j else 15 FIXED AREA SIZE so set up one boundary xMax size float 2 xMin gt xMax yMin xMin yMax xMax zMin xMin zMax xMax if inCFile bad cerr lt lt Unable to read in all configuration values lt lt endl lt lt Some variables may use default values or even cerr unsuitable values endl inCFile close auxInitDisplayMode AUX SINGLE AUX Set initial window position on screen auxInitPosition lowerLeft X lowerLeftY upperRightX upperRightY viewport 0 viewport 1 viewport 2 viewport 3 auxInitWindow argv 0 Opens a window
46. and top of the cylinder Then the cylinder is lengthened and trimmed using a cut plane calculated from the two points that surround the point used for each end see Figure 1 This cut plane is defined as the plane perpendicular to the line connecting the E two surrounding points and containing the end point The OpenGL procedure gluCylinder is used to build each cylinder In OpenGL cylinders can have different radii at each end However since gluCylinder creates cylinders with one end at the origin and the axis along the positive z axis each piece of the loop must first be translated so that its initial end is at the origin and then rotated so that its axis lies along the z axis Then the cylinder can be created The cylinder is lengthened as mentioned above to make sure that the cut planes lie entirely within the cylinder is center of end plane at left ofcone B is center of end plane at right of cone C is center of end plane for cone to right of this cone L is center line of cone P is perpendicular to line from to C R is radius of plane at B Plane at B is defined by B lies on plane and it is perpendicular to Figure 1 Tapered Cylinder Building Block After the cylinder is created it is placed back in the image where it originally was by reversing the rotation and translation operations This reversal is planned for by saving the projection matrix before moving the cylinder to t
47. apFileHeader bitMapFileHeader bfType 0x4D42 bytes bitMapFileHeader bmiHeader colortable image bitMapFileHeader bfSize 14 40 256 4 WIDTH HEIGHT bitMapFileHeader bfReservedl 0 bitMapFileHeader bfReserved2 0 BYTES bmiHeader colortable bitMapFileHeader bfOffBytes 14 10 4 256 4 struct 1 DWORD biSize of bytes required by this structure DWORD biWidth width of bitmap in pixels DWORD biHeight height of bitmap in pixels WORD biPlanes of planes for target must 1 WORD biBitCount of bits per pixel DWORD biCompression Type of compression 0 for no compression DWORD biSizelmage of bytes in the image DWORD biXPelsPerMeter Horizontal resolution in pixels per meter DWORD biYPelsPerMeter Vertical resolution in pixels per meter DWORD biClrUsed color table reference 0 for immediate table DWORD biClrImportant set to O if all colors are important bmiHeader bmiHeader biSize 40 bmiHeader biWidth WIDTH bmiHeader biHeight HEIGHT bmiHeader biPlanes 1 bmiHeader biBitCount 8 bmiHeader biCompression 0 bmiHeader biSizelmage WIDTH HEIGHT bmiHeader biXPelsPerMeter 3000 bmiHeader biYPelsPerMeter 3000 bmiHeader biClrUsed 0 bmiHeader biClrImportant 0 imageOut setf ios binary Put the bitMapFileHeader in the bmp file memopy buff amp bitMapFileHeader bfType 2 2 amp bitMapFileHeader
48. center xVal top center xVal cylVect yParam bot center y Val top center y Val cylVect zParam bot center zVal top center z Val cylVect offset 0 0 Define the bottom cut plane CutPlane prior center bot center top center amp cutPl1 dot DotProd cylVect cutPI1 angZ1 RadToDeg float acos fabs dot CheckAngle angZ1 bot center Make sure angle is within limits Calculate cut plane for top CutPlane bot center top center after center amp cutPI2 dot DotProd cylVect cutPI2 angZ2 RadToDeg float acos fabs dot CheckAngle angZ2 top center Make sure angle is within limits SetEqn cutP11 eqnl 1 glClipPlane GL CLIP PLANEO eqnl clip for bottom end of cylinder glEnable GL CLIP PLANEO SetEqn cutPI2 2 1 glClipPlane GL CLIP PLANEI eqn2 clip for top end of cylinder Clip planes are in eye coordinates not scene coordinates so they will be a rendered a little bit off what they really are The mathematical equations will be correct and the resulting x ray image will be correctly calculated 58 c T c V AU IEEE glEnable GL CLIP PLANE Hmc Position the cylinder Translate so that bottom end is at origin Rotate so that centerline is on Z axis Calculate bottom cut plane Raise above Z axis by enough that the angled cut plane will not go below XY plane before it finishes cutting the end Af
49. cylinder It also takes two more DISKs containing the prior and after points and radii It makes a tapered cylinder using the DISKs for bottom and top of the cylinder Then the cylinder is lengthened and trimmed using a cut plane calculated from the two points that surround the point used for each end This cut plane is defined as the plane perpendicular to the line connecting the two surrounding points and containing the end point MakeCylinder also calls MapCylinderPixels to place depth info about this cylinder in the pixel array This needs to be done at the time that the complete projection matrix for this cylinder is available void MakeCylinder DISK prior DISK bot DISK top DISK after PLANE cutPl1 cutPI2 cylVect GLfloat angZ1 angZ2 deltaZBot 0 0f 0 0f dot DISK newTop GLdouble eqn1 4 GL double eqn2 4 GLdouble density glColor3f 0 66 0 06 0 40 shade of violet glPushMatrix Save matrix position jy ecco a cM HAS UNE m TONES X Calculate the angle between the cut planes and the axis of the cylinder If the plane is perpendicular to the axis of the cylinder angZ is 0 deg cos angle u v u v where and v are vectors plane normal v cylinder axis Thus angle arccos u v u v in radians Use the absolute value of the dot product so that angle is always positive Cylinder vector centerline cylVect xParam bot
50. e is drawn The three EYE values are the coordinates of the viewing point The three CENTER values are the coordinates of the center of the scene The three UP values determine the direction shown as the top of the scene on the screen The CENTER values are designed to be the origin The EYE values work best if they are each less than 1 The CENTER and EYE values really only determine the direction from which the scene is viewed See the explanation for line 4 below Line 3 controls the x ray picture generation The sensitivity function index M the sensitivity function to be used This relates to the type of filter in use The exposure value sets the exposure which relates to the brightness of the image Pixel sizing defines the type of scaling to be used FIXED AREA SIZE will show only that part of the image that is within the fixed boundaries FIXED PIXEL SIZE will show the image as it would look with each pixel covering the stated amount of area on the Sun s surface If line 3 has FIXED AREA SIZE then size contains the horizontal span for the fixed area the y and z values match the x one to make the viewing space a cube If line 3 has FIXED PIXEL SIZE then the value of size is the span of one pixel Line 4 is the offset 18 for the center of the viewing volume The x y and z values entered here will be the center of the volume NOTE length parameters radius x y Z pixel size and offset must be in the same units Positive z
51. element in column 1 in the top position Next divide the top row by the first element This makes the first element 1 The top row ofthe identity matrix is also divided by this element Now all other rows in the matrix are individually added to a multiple of the first row chosen so that the resulting first element will be zero After completing this operation the entire first column has been reduced to a 1 on the diagonal the first element in this case and zeroes for all other elements This sequence is now repeated for each of the rest of the columns The transformation matrices in openGL are stored in column major order This means that A 1 is the first element of the second row A 4 is the second element of the first row ie A 0 A 4 12 5 A 9 13 A 2 A 6 A 10 A 14 A 3 A 7 11 A 15 FECES CT ICICI IIIA void InvertMatrix GLdouble inputMatrix DIM DIM GLdouble outputMatrix DIM DIM int row col k int pivot GLdouble mat DIM DIM outMat DIM DIM GLdouble factor Copy inputMatrix to temp storage for k 0 k lt DIM DIM k mat k inputMatrix k Initialize outMat to the identity matrix for 0 k lt DIM DIM k Make all zero outMat k 0 for row 0 row lt DIM row Set diagonals to 1 71 outMat row DIM row 1 Begin G J operations for row 0 row lt DIM rowt Walk down all rows
52. erature value into a sensitivity value using the sensitivity function specified New functions will be added to this switch statement as they are needed void ApplySensitivityFcn int sensFcn GLdouble temp switch sensFcn case 0 if temp lt 1 0e 6 temp 0 0 else temp 1 0 break case 1 60 temp lt yohkoh instr rsp 00 temp break case 2 temp yohkoh instr rsp 01 temp break case 3 temp yohkoh instr rsp 02 temp break case 4 temp yohkoh instr rsp 03 temp break case 5 temp yohkoh instr rsp 04 temp break case 6 temp lt yohkoh instr rsp 05 temp break default if temp lt 1 06 temp 0 0 else temp 1 0 Routines for T response of Yohkoh given T in degrees K returns f T in 10 28 DN sec YkP BCD Yohkoh pixel pu BCD brems col depth int n e 2 dl in cm 5 routines include yohkoh instr rsp 00 open yohkoh instr rsp 01 Al 1 yohkoh instr rsp 02 Al Mg yohkoh instr rsp 03 Mg3 yohkoh instr rsp 04 2 yohkoh instr rsp 05 86119 static double yohkoh instr rsp 00 double temp response code for Yohkoh open int i double lgt frc static double dlog10t 0 0500002 static double tv 51 5 500006 00 5 55000e 00 5 60000e 00 5 65000e 00 5 70000e 00 5 75000 00 5 80000e 00 5 85000 00 5 90000 00 5 95000e 00 6 00000 0
53. ffer pixel color proportional to the distance of the coordinate transformations are accomplished by matrix multiplication The transformation from cylinder coordinates to screen coordinates uses just the projection transformation already calculated by OpenGL The inverse of that matrix is used to go from screen to cylinder coordinates The inverse matrix is calculated by the Gauss Jordan Elimination method Press 32 37 The surface of a tapered cylinder can be represented by a mathematical construct known as a quadric This isa 4 x 4 matrix The values in the matrix are determined from the equation of the surface X x y y r1 r1 2 1 1 2 2 2 2 r1 r2 z zr2y 0 13 where r1 and r2 are the radii at the two ends of the tapered cylinder The cylinder has one end at 2 0 and the other end at z 272 The general form of the quadric equation is Kirk 280 F X y z 2bxy 2cxz 2dx ey 2fxy 2gy hz 212 0 The matrix coefficients are Q11 012 021 b 013 Q31 014 041 d 022 e etc Fora right tapered cylinder with z axis of symmetry values of the non zero elements are Q11 022 1 033 r1 2 2 2 Q34 043 rl r112y h and Q44 r1 Where r1 is the radius at 2 0 r2 is the radius at z h and h is the height of the cylinder So only Q33 Q34 and Q44 need to be calculated for each cylinder To calculate the intersections we need
54. for Yohkoh Bel 19 int i double lgt frc static double 010810 0 0500002 static double tv 51 5 50000 00 5 55000e 00 5 60000e 00 5 65000 00 5 70000e 00 5 75000e 00 5 80000e 00 5 85000e 00 5 90000e 00 5 95000e 00 6 00000 00 6 05000 00 6 10000e 00 6 15000 00 6 20000e 00 6 25000 00 6 30000 00 6 35000 00 6 40000 00 6 45000 00 6 50000 00 6 55000 00 6 60000 00 6 65000 00 6 70000 00 6 75000 00 6 80000 00 6 85000 00 6 90000 00 6 95000 00 7 00000 00 7 05000 00 7 10000 00 7 15000 00 7 20000 00 7 25000 00 7 30000 00 7 35000 00 7 40000 00 7 45000 00 7 50000 00 7 55000 00 7 60000 00 7 65000 00 7 70000e 00 7 75000 00 7 80000 00 7 85000 00 7 90000 00 7 95000 00 65 8 00000e 00 static double fv 51 1 62282e 18 8 69929e 17 3 51845e 15 1 03648e 13 2 41363e 12 4 21862e 11 1 63456e 09 1 54806e 08 1 48222e 07 1 06793e 06 7 02349e 06 3 46967e 05 1 41041e 04 4 99637e 04 1 62981e 03 4 89867e 03 1 32784e 02 3 15088e 02 6 68034e 02 1 30267e 01 2 05412e 01 3 34308e 01 5 20663e 01 7 78464 01 1 12565e 00 1 57786 00 2 14043 00 2 79409 00 3 52412 00 4 34205 00 5 15012 00 5 90210 00 6 55741 00 7 07878 00 7 49269 00 7 82844 00 8 12285 00 8 39323 00 8 65198 00 8 90404 00 9 15081 00 9 38752 00 9 60492 00 9 79802 00 9 96492 00 1 01053 01 1 0218
55. g to be used FIXED AREA SIZE will show only that part of the image that is within the fixed boundaries whereas FIXED PIXEL SIZE will show the image as it would look when each pixel covered the stated amount of area on the Sun s surface Ifline 3 has FIXED AREA SIZE then size contains the horizontal span for the fixed area size the y and z values match the x ones and thus all viewing is of a cube If line 3 has FIXED PIXEL SIZE then the value is the span of one pixel Line 4 is used to offset the center of the viewing volume Place the coordinates of the desired center for the viewing volume here Be sure to use the same units as for all other dimensions NOTE All length parameters radius x y z pixel size and offset must be in the same units Positive Z is away from the Sun s surface File input data 4 58318 5 03010 2 39816 1 39698 32 5637 318820 4 58318 5 75791 2 96428 0 00000 32 5637 318820 4 35026 5 03010 2 39816 1 39608 32 5637 928480 4 23094 4 31004 1 89874 2 35220 13 6805 1 43248 06 4 12981 3 58737 1 36782 3 28825 8 03442 1 86923 06 4 04729 2 84821 0 804874 4 19207 5 53772 2 25151 06 3 98117 2 08341 0 211612 5 05411 4 18923 2 58864 06 3 92774 1 29056 0 407562 5 87134 3 36101 2 89004 06 3 88476 3 85092 3 82305 3 79963 3 77955 3 76126 3 74373 3 72613 3 70771 3 68810 3 67124 3 66537 3 64236 3 61783 3 59181 3 56476 3 53735 3 51022 3 48440 3 46086 3 44073 3 42629 3 41876 3 41983
56. gree that the Library shall make it available to borrowers under rules of the Library If I have indicated my intention to copyright this thesis by including a copyright notice page copying is allowable only for scholarly purposes consistent with fair use as prescribed in the U S Copyright Law Requests for permission for extended quotation from or reproduction of this thesis in whole or in parts may be granted only by the copyright holder Signature _ Date iv TABLE OF CONTENTS TM vii 1 Problem Background i oes uoa aset eese mh ee a e on dep Solar Magnetic Phenomena 2 Understanding Solar Processes NOUS cd 3 SXI REQUIREMENTS sese MET 4 SXIINTERNALS PTT Building the Plasma Loop Mapping to Pixel 00105 ROREM 8 Mapping Transformations S 9 9 14 CONCLUSION AND POSSIBLE IMPROVEMENTS 15 16 REFERENCES CITED 19 APPENDICES 20 APPENDIX A SALCPP UPPER du 217 APPENDIX B CYLINDE
57. he origin After creating the cylinder simply restoring the previous projection matrix puts the cylinder back where it belongs For each set of four records in the input data file one cylinder can be drawn The process proceeds by sliding over one record and producing another cylinder until the end of the data file is reached This means that each record in the data file except the ones at the beginning and end is used as an end twice for each of the cylinders that meet at that point and four times to help define a cut plane Mapping to Pixel Colors After each tapered cylinder is drawn it is mapped to pixel colors The color of each pixel is proportional to the depth density and temperature of the portion of the plasma loop that maps to that pixel Here depth means the distance through the plasma loop as seen from the pixel in question The process requires several steps To begin with this is an area where designing for speed is critical We are about to map this cylinder on a pixel by pixel basis and do not want to deal with pixels that do not map from this cylinder The depth is the distance through the cylinder as experienced by a ray emanating from the pixel on the screen and traveling perpendicular to the screen The density value 9 used is actually the square of the electron number density divided by 10 which is read directly from the input data file The temperature value is a function of the plasma temperature
58. hen update values of min and max in screen coords Loop over all eight points ILL DEMO ant nt glGetDoublev GL MODELVIEW MATRIX modelMat glGetDoublev GL PROJECTION MATRIX projMat MatrixMult initInverse modelMat modelMatFixed init TRUE Used to initialize minMax Bottom points 2 0 cylPt zVal GLfloat 0 0 glColor3f 0 0f 0 0f 1 0 blue for i 1 i lt 2 2 cylPt xVal i SQRT2 botRad for j 1 j lt 2 2 42 cylPt yVal 1 SORT2 botRad gluProject cylPt xVal cylPt y Val cyIPt zVal modelMat projMat viewport amp scrPt xVal amp scrPt yVal amp scrPt zVal MinMax minMax scrPt amp init Update x and y min max values Top points z height cylPt zVal height for i 1 i lt 2 2 cylPt xVal 1 SQRT2 topRad for j 1 j lt 2 j 2 cylPt yVal j SORT2 topRad gluProject cylPt xVal cylPt yVal cyIPt zVal modelMat projMat viewport amp scrPt xVal amp scrPt yVal amp scrPt zVal MinMax minMax scrPt amp init Update x and y min max values Part 3 Part 4 and Part 5 included in this section Now minMax contains points representing the bounding rect in screen pixel coordinates Part 6 Set the color for the screen pixels This will be done using an additive process so that if other cylinder points are mapped to this pixel their contribution to the total can be
59. is away from the Sun s surface Examples of sxi ini and a data file are included in Appendix G An example of the program command line would be sxi input_data where input_data is the name of the file containing the loop data and the file sxi ini is in the curent working directory REFERENCES CITED Foley VanDam Feiner Hughes and Phillips Introduction to Computer Graphics Reading Massachusetts Addison Wesley Publishing Company 1995 Kirk D ed Graphics Gems III Boston Harcourt Brace Jovanovich c1992 Lockheed First Light online Available http www space lockheed com SXT html2 First_Light html March 24 1998 Longcope D W Topology and Current Ribbons A Model for Current Reconnection and Flaring in a Complex Evolving Corona Solar Physics 169 1996 91 121 Press W H Numerical Recipes in C The Art of Scientific Computing Cambridge Cambridge University Press 1988 Reale and Peres G Solar Flare X ray Imaging Coronal Loop Hydrodynamics and Diagnostics of the Rising Phase Astronomy and Astrophysics 299 1995 225 237 Woo Neider and Davis OpenGL Programming Guide Second Edition Reading Massachusetts Addison Wesley Developers Press c1997 APPENDICES 21 APPENDIX A SXI CPP 22 sxi cpp Steve Lundberg 2 15 97 This is the driver file for viewing the intensity of x rays produced by solar flares defined by solar m
60. k void MapCylinderPixels GLdouble botRad GLdouble topRad GLdouble height PLANE PLANE cutPI2 GLdouble density GLdouble sensitivity1 GLdouble sensitivity2 41 int init PLANE scrPlane POINT3 cylPt scrPt int i j GLdouble modelMat DIM DIM GLdouble projMat DIM DIM GLdouble modelMatFixed DIM DIM GLdouble dist int counts GLdouble minMax 4 index 0 is x min 1 is x max 2 is y min 3 is y max int xPos yPos 1 vector eye pt center pt scrPlane xParam GLfloat EYEX CENTERX scrPlane yParam GLfloat EYEY CENTERY scrPlane zParam GLfloat EYEZ CENTERZ Part 2 The bounding rect for the cyl is easily found in the coords where it is initially created Using the un clipped cyl gives slightly larger values than absolutely necessary but is much easier When the cyl is created its bottom lies in the x y plane z 0 and extends out a radius When viewed from any angle the points x y SORT2 radius will always be bounds for this circle Similarly for the top So if the eight points 4 for the bottom and 4 for the top are transformed to screen coords and then used to generate min max values the resulting rect bounds the cylinder on the screen Using gluProject will give the x and y values for this rectangle in pixel numbers from the origin of the screen lower left corner Calculate the transform of one of the eight points t
61. le inFile eof amp amp dataCheck END PLASMA LOOP CopyDisk testDisk1 amp testDiskO CopyDisk testDisk2 amp testDisk1 CopyDisk testDisk3 amp testDisk2 inFile getline dataLine MAX DATA LINE LENGTH if dataCheck CheckDataLine GOOD DATA dataLine amp testDisk3 radius amp testDisk3 center xVal amp testDisk3 center y Val amp testDisk3 center z Val amp testDisk3 density amp testDisk3 temperature Apply sensitivity function to temperatures of testDisk1 and testDisk2 This will give the value that will actually be used in deciding the color of the pixels from this sylinder ApplySensitivityFcn sensFcn amp testDisk1 temperature ApplySensitivityFcn sensFcn amp testDisk2 temperature MakeCylinder testDisk1 testDisk2 testDisk3 of loop to make all cylinders in one plasma loop inFile eof dataCheck END OF FILE return dataCheck Create an image file in bmp format FECA CISC COIS IE A CICS A IAC TTT void Savelmage GLubyte xrayData HEIGHT WIDTH ofstream imageOut char buff 80 char color 4 58 struct WORD bfType file type must be BM DWORD 5125 size of this file in bytes WORD bfReservedl reserved and must 0 WORD bfReserved2 reserved and must 0 DWORD bfOffBytes byte offset from start to the actual bitmap bitM
62. matical model of solar dynamics The time it takes is short enough to support making movies of the time development of these plasma loops The input is given in units typically used for solar modeling It supports plasma loops where the temperature and density vary along the axis of the loop The pixel size can be fixed at some desirable size or allowed to adjust to make the plasma loop fit in the viewing area It synthesizes images using different camera filters and sensitivities It is portable between Windows and Unix environments Several modifications may enhance the usability of SXI o Currently the viewing volume is centered on the origin Removing this requirement would allow easier interface to model output Some textual output on the synthetic image would be desirable This might include data identification filter in use and The coordinates could be changed from to spherical to be more consistent with the usual solar coordinate use SXI USER S MANUAL SXI uses two input files the data file and the initialization file The data file must be named on the command line while the initialization file has a fixed name and must reside in the current working directory The data file format is one six entry record per line The six numeric values in each record are radius x coordinate y coordinate z coordinate density value and temperature The radius and coordinate
63. matrix so can map from screen coordinates to current cylinder coordinates First retrieve the matrix then invert it The inverted matrix is returned at the supplied pointer From Press W H Numerical Recipes in C Second Edition Gauss Jordan Elimination G J Elim is as good as any other technique for inverting matrices For the closely related case of solving sets of linear equations LU decomposition is better In any matrix inverting operation one needs to be careful of singular matrices Initially I will not support singular matrices but the technique of Singular Value Decomposition SVD discussed in Press will handle some kinds of singular matrices G J Elim with pivoting starts with the upper left element of the matrix 70 An equally sized identity matrix is also used Operations on the initial matrix are also performed on this identity matrix which is then no longer the identity matrix The goal of this procedure is to transform the initial matrix into an identity matrix by using only linear operations When the initial matrix has been transformed to the identity matrix the original identity matrix is the inverse matrix of the initial matrix The pivoting is performed to guarantee the mathematical stability of this operation If pivoting is not used there is a substantial probability of roundoff and truncation errors dominating the solution obtained First pivot the rows of the matrix to place the largest
64. nleft Angle bisector is the plane perpendicular to line from left to right only if the length from the middle to each of the other points is the same So first find a new right point that is on original vector to right point but is exactly as far away as the left point is from the middle Beh RTT NE Y lenright float sqrt pow right xVal middle xVal 2 pow right yVal middle yVal 2 pow right zVal middle zVal 2 lenleft float sqrt pow left xVal middle xVal 2 pow left yVal middle yVal 2 pow left zVal middle zVal 2 Now construct vector from new right to original left points This is The normal to a plane that is perpendicular to the vector The equation for the vector from new right to original left point is D A Where D is end of vector from B in direction of original C B vector and of same length as vector from A to B or D B new right vector The new right vector starts at B goes in the direction C B and is adjusted to have same length as vector from A to B B A C B Thus D is B This yields the equation 36 D A B A Where is left point B is middle point is original right point D is new right point B A is length of left vector and C B is length of original right vector This equation holds for each of the three directions pl gt xParam middle xVal
65. nother requirement of the program is that it run on as many platforms as possible so that it may be used by many researchers The images generated are close approximations to the actual plasma loops described by the models They also look much like actual x ray images taken of the solar plasmaloops The program is able to generate an image of several plasma loops in less than one minute This allows the generation of an entire movie consisting of perhaps a hundred images in a few hours Because C and OpenGL are available on both Windows and Unix platforms the program can be compiled and run on many of the machines used by researchers around the world The program works as intended and will be a substantial help in the further development of the solar magnetic models It is general enough that it can be used by many different researchers using different solar magnetic models INTRODUCTION Problem Background Physicists are studying the Sun in an attempt to understand the processes at work there Longcope One of the tools they use is x ray pictures taken of the corona The plasma in the corona is hot enough that it emits x rays so features in the corona such as the plasma loops there can be seen in x ray photographs The plasma loops in the corona are driven by convection and other transport processes occurring within the Sun Some understanding of the processes at work below the surface of the Sun can be gained by studying the loop
66. om the interior of the Sun Often these features are connected by magnetic field lines within the corona They influence each other through these connections One possible configuration is with two magnetic features opposing each other e g North to North The features can be driven toward each other by the underlying dynamics within the Sun but they repel each other through their magnetic field interaction Energy can be stored in this interaction When the energy stored reaches some threshold the two magnetic fields reconnect their field lines in a new combined configuration This releases much of the energy stored in the previous configuration This energy gets deposited as heat in the plasma where the field lines are The increased heat in the plasma makes it emit more x rays These x ray emissions are the probe being used to follow the magnetic field interactions Plasma is in the magnetic field lines above the Sun s surface because as the magnetic field lines are pushed up 3 from within the Sun they drag the plasma along with them This occurs because charged particles cannot readily cross magnetic field lines Understanding Solar Processes One common way to formulate our understanding of physical processes is to build models of those processes When these models closely approximate what happens in nature we believe the models to be accurate In the case of understanding solar dynamics models of the formation and time devel
67. opment of plasma loops and systems of interacting plasma loops have been developed However the models are numerical models and it is difficult to compare the tables of model output with x ray images taken of plasma loops on the Sun One very powerful method for comparing the models to actual events is to create synthetic x ray pictures of the loops in the models and compare these pictures to the ones taken of the Sun However at present only relatively crude methods are used to generate pictures of these model plasma loops Reale The value of good color renditions of the plasma loops generated by the models is readily apparent The opportunity to develop a program that would be of great value to Dr Acton s group and others in the solar physics community resulted in the Synthetic X ray Imager SXI computer program described in this thesis SXI REQUIREMENTS The primary goal of the SXI program is to support the development of solar magnetic models As such it needs to have an easy and natural interface SXI is designed to use the same parameters and measurement units as the solar magnetic models Since SXI is intended to create images that look similar to the ones taken by the X ray camera on the Yohkoh satellite it also uses parameters and units typical of those used on that satellite The solar magnetic models describe the plasma loops they are dealing with by 8 series of points and associated attributes Each point in space
68. put cylinder back where it belongs glBegin GL LINES Draw centerline of cylinder xem glColor3f 0 0 1 0 1 05 glVertex3d bot center xVal bot center y Val bot center z Val glVertex3d top center xVal top center y Val top center z Val glEnd Move this point of the object to the origin void MoveToOrigin POINT3 pt glTranslated pt xVal pt y Val pt zVal Calculate rotation vector and rotation angle to bring centerline to Z axis void AlignZAxis POINT3 20113 topPt 35 POINT pt GL float ang Set pt end of vector from origin parallel to centerline pt xVal topPt xVal botPt xVal pt yVal topPt yVal botPt yVal pt zVal topPt zVal botPt zVal Calculate degrees to rotate angle between vector from origin to pt and the Z axis Nength sqroot of pt xVal 2 pt yVal 2 pt zVal 2 angle lt arccos z length ang GLfloat RadToDeg acos pt zVal sqrt pow pt xVal 2 pow pt yVal 2 pt zVal 2 glRotatef ang GLfloat pt y Val GLfloat pt xVal GLfloat 0 0 RAR HR RAE RAH A RR HH HHH HH HH HH HH HHH Return plane through middle point and bisecting the angle made by the lines from the left and right points Planar form Ax By Cz D 0 void CutPlane POINT3 left POINT3 middle POINT3 right PLANE pl float lenright le
69. re is no intersection void IntersectRayWithCutPlane POINT3 rootPtr PLANE cutPl RAY scrRay GLdouble t denom cutPl xParam scrRay dirVect x Val cutPl yParam scrRay dirVect yVal cutPl zParam scrRay dirVect zVal if denom 0 0 Ray is parallel to plane and outside of cut cylinder This should never happen since we have already checked printf Ray did not intersect cut plane exit IMAGINARY INTERSECTION 51 else t cutPl xParam scrRay initPt xVal cutPl yParam scrRay initPt yVal cutPl zParam scrRay initPt zVal cutPl offset denom rootPtr gt xVal GLfloat scrRay initPt xVal scrRay dirVect xVal t rootPtr gt GLfloat scrRay initPt y Val scrRay dirVect yVal t rootPtr gt zVal GLfloat scrRay initPt zVal scrRay dirVect zVal t Set color of designated pixel proportional to distance through cylinder Dist is scaled before it is sent here The incoming x and y positions are screen pixel numbers they need to be scaled to fit the pixel buffer which may not be the full screen size void AdjustPixelColor int xPos int yPos int dist GLubyte colorValue Negative distance shouldn t happen if dist lt 0 dist 0 cout Ray distance thru cylinder was negative x y xPos yPos endl xPos xPo
70. rtho xOrthoMin xOrthoMax yOrthoMin yOrthoMax zOrthoMin zOrthoMax E MEHR ROTEN TORT 1 Origin is at middle of screen if xMin xMax and yMin yMax Ortho view mode makes cut cube from the 6 planes defined in x y z min max order Viewing is from the negative z axis This can be modified by using the gluLookAt function Because this is a parallel projection system the distance from the viewer to the viewing volume is immaterial All that matters is mapping the volume to screen coordinates and setting the direction of view Ae a ELE IO Eric SNC glMatrixMode GL MODELVIEW glLoadIdentity glClearColor 0 52 0 5f 0 5f 0 07 set background color glClear GL COLOR BUFFER DEPTH BUFFER gluLookAt EYEX EYEY EYEZ CENTERX CENTERY CENTERZ UPX UPY UPZ Perhaps rather than making a string of cylinders in which the matching of the ends is rather difficult and unavoidably contains some error it would be possible to use a nurbs surface created from the necklace of points This surface would probably need some tweaking to get it close enough to the points but it might work quite well It certainly would eliminate cylinder end point problems This idea will not be pursued for now FEAE k aae ak ae a ae ake ae ake ae ake ae ake ae ae ae ake ae a ake ae ae ae ake ake ae ae a ae ake
71. s WIDTH viewport 2 yPos yPos HEIGHT viewport 3 if lt HEIGHT amp amp xPos lt WIDTH yPos gt 0 amp amp xPos gt 0 accumulate into pixel image buffer dist int xrayData yPos xPos Saturation is 255 if dist gt 255 dist 255 colorValue GLubyte dist xrayData xPos colorValue else cout lt lt Bad xrayData index y x lt lt yPos lt lt lt lt xPos lt lt endl 52 APPENDIX D SXIUTILS CPP sxiUtils cpp Steve Lundberg 2 15 97 Defines the function used to draw the synthetic xray picture Defines the functions used to read in the data file check the data for validity and start the cylinder creation process The function to output the bitmap file is here For future possible use Uses the rand fen to get numbers between 0 and 255 These values can be used to produce a surface plane to represent the surface of the sun include lt stdlib h gt include lt fstream h gt include lt process h gt include lt math h gt ifdef WIN32 include lt afxwin h gt include lt wingdi h gt endif include lt GL gl h gt include lt GL glu h gt include lt GL glaux h gt include sxi h Local prototypes static double yohkoh instr rsp 00 double temp static double yohkoh instr rsp 01 double temp static double yohkoh instr r
72. s is for the near clipping plane Find the start point of this ray in cylinder coordinates gluUnProject scrPt xVal scrPt yVal scrPt zVal modelMat projMat viewport amp cylStartPt xVal amp cylStartPt yVal amp cylStartPt z Val scrPt zVal 2 0 This is for the far clipping plane Find the ending point of this ray in cylinder coordinates gluUnProject scrPt xVal scrPt y Val scrPt zVal modelMat projMat viewport amp cylEndPt xVal amp cylEndPt y Val amp cylEndPt zVal Calculate the vector from the startPt to the endPt in cylinder coordinates cylVect xVal cylEndPt xVal cylStartPt xVal cylVect yVal cylEndPt y Val cylStartPt y Val cylVect zVal cylEndPt zVal cylStartPt zVal Make the direction vector have unit length lenVect sqrt pow cylVect xVal 2 pow cylVect yVal 2 pow cylVect zVal 2 cylVect xVal cylVect xVal lenVect cylVect yVal cylVect yVal lenVect cylVect zVal cylVect zVal lenVect cylRay initPt cylStartPt cylRay dirVect cylVect jJ c MEIN Lun c OM m Solve equation for roots intersections Need equation for the cylinder in question in cylinder coordinates this is the quadric The equation is in the form of k2 t t k1 t 0 For the case of a circular truncated cone called cylinder for short k2 ax ax Q11 ay ay Q22 az az Q33 2 ax x0 Q11 2 ay y0 Q22 2 az z0 Q33 Q34 x0 x0 Q11 y0 y0 Q2
73. s must be given in the same units but the choice of units is arbitrary The density sinu the square of the electron number density divided by 10 which is expected to be a number between zero and 100 The temperature is given in degrees Kelvin The values must be separated by white space and terminated with an end of line They will be read in floating point format and can be in decimal notation or exponential base ten notation If more than one plasma loop is contained in file the special string LOOP must be placed on the line immediately following the last record for each loop The next record will be taken as the first record for a new loop The first and last records for each loop are only used to determine the cut plane for the ends of the plasma loop They do not contribute to the loop itself The initialization file is named sxi ini It has four lines of setup information The file format is line 1 windowX windowY windowWidth windowHeight 17 line 2 EYEX EYEY EYEZ CENTERX CENTERY CENTERZ UPX UPY UPZ line3 Sensitivity function index exposure value pixel sizing 5126 line 4 xOffset yOffset zOffset Line 1 controls the position and size in pixels of the window on the screen This needs to be the same size as the final synthetic x ray image or pixel artifact effects will be seen Currently the synthetic x ray image is 600 by 600 pixels Line 2 controls the viewpoint from which the scen
74. s valid process it as part of the current plasma loop If the line is END OF LOOP code then prepare to render another plasma loop At end of file save resulting picture in a suitable data file preferably gif format FUTURE An optional title may also be entered in the data file The 30 title is on the line that begins with the code TITLE CEU dataCheck GOOD DATA while dataCheck END OF FILE dataCheck MakePlasmaLoop inFile inFile close glPopMatrix end of make plasma loops section Display the xray image pixel buffer statusOK displayData xrayData Save the synthetic xray image as a BMP file Savelmage xrayData imageOut glFlush 31 APPENDIX B CYLINDER CPP 32 cylinder cpp Steve Lundberg 10 10 97 This file contains functions for drawing a tapered cylinder include lt fstream h gt include lt process h gt include lt stdlib h gt include lt math h gt include lt stdio h gt ifdef WIN32 include lt afxwin h gt include lt wingdi h gt endif include lt GL gl h gt include lt GL glu h gt include lt GL glaux h gt include sxi h GLOBAL VARIABLES GLUquadricObj cylQuadric Used for the construction of each cylinder MakeCylinder takes two DISKs containing the center position and radius at the bottom and top of the
75. s visible in an x ray picture Since x rays do not penetrate the Earth s atmosphere the x ray pictures of the Sun must be taken from satellites The work described here was done in conjunction with the MSU Solar Physics group led by Dr Loren Acton Dr Acton s group is one of several research organizations using the Yohkoh satellite to take x ray pictures of the Sun The Yohkoh satellite is the product of a joint Japanese U S venture This satellite is operated by Lockheed and current generally icis 68 24 hours old x ray pictures of the Sun are posted on the Internet Lockheed Dr Acton then at Lockheed directed the design and construction of the Soft X ray Telescope SXT carried on Solar Magnetic Phenomena Much of what can be seen in x ray pictures of the Sun is believed to be the result of magnetic fields interacting with the plasma in the Sun s corona These magnetic fields are believed to be driven by electrical currents within the Sun The electrical currents are the result of convective processes occurring inside the Sun As the magnetic fields interact with the coronal matter the matter is transported and heated The pattern of heating determines the quantity and distribution of x ray emission Magnetic field enters the corona from the interior of the Sun through isolated magnetic features on the solar surface These features are the top part of magnetic features called flux tubes fr
76. sp 02 double temp static double yohkoh instr rsp 03 double temp static double yohkoh instr rsp 04 double temp static double yohkoh instr rsp 05 double temp A69 che e ke oe ke ke ke ke ok ok oe je e ke he oj je ne ke e oj he je je ke ke ok Draws the pixel data in the array at the current screen location ATA AR RAAT A AHA AAA int displayData GLubyte xrayData HEIGHT WIDTH glDrawPixels WIDTH HEIGHT GL COLOR INDEX GL UNSIGNED BYTE xrayData return 0 void SetupAxes void 54 glPushMatrix glBegin GL LINES Draw the three axes glColor3f 1 0f 0 0f 0 0f red glVertex3d xMin 0 0 glVertex3d xMax 0 0 glColor3f 0 0f 1 0 0 0f green glVertex3d 0 yMin 0 glVertex3d 0 yMax 0 glColor3f 0 0f 0 0f 1 0f blue glVertex3d 0 0 zMin glVertex3d 0 0 zMax glEnd glPolygonMode GL FRONT AND BACK GL LINE Draw 3 squares in the x y plane 886 GL POLYGON glColor3f 0 0f 0 0f 1 0f blue glVertex3d xMin 3 yMin 3 0 glVertex3d xMin 3 yMax 3 0 glVertex3d xMax 3 yMax 3 0 glVertex3d xMax 3 yMin 3 0 glEnd glBegin GL POLYGON glColor3f 0 0f 0 00 1 0f blue glVertex3d 2 xMin 3 2 yMin 3 0 glVertex3d 2 xMin 3 2 yMax 3 0 glVertex3d 2 xMax 3 2 yMax 3 0 glVertex3d 2 xMax 3 2 yMin 3 0 glEnd glBegin GL POLYGON glColor3f
77. sting of perhaps a hundred images in a few hours Because C and OpenGL are available on both Windows and Unix platforms the program can be compiled and run on many of the machines used by researchers around the world The program works as intended and will be a substantial help in the further development of the solar magnetic models It is general enough that it can be used by many different researchers using different solar magnetic models SYNTHETIC X RAY IMAGER FOR SOLAR PLASMA LOOPS by Steven Kenneth Lundberg A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science MONTANA STATE UNIVERSITY Bozeman Montana May 1998 191 APPROVAL of a thesis submitted by Steven Kenneth Lundberg This thesis has been read by each member of the thesis committee and has been found to be satisfactory regarding content English usage format citations bibliographic style and consistency and is ready for submission to the College of Graduate Studies 4 Date Chairpersofi Graduate Committee Approved for the Major Department 1 7 1 0 Data Head Major Department Approved for the College of Graduate Studies FE Date Graduate D 111 STATEMENT OF PERMISSION TO USE In presenting this thesis in partial fulfillment of the requirements for a master s degree at Montana State University I a
78. t x Val posRootPtr gt cylRay initPt yVal GLfloat posRootT cyIRay dirVect y Val posRootPtr gt zVal cylRay initPt zVal GLfloat posRootT cylRay dirVect zVal negRootPtr gt xVal cylRay initPt xVal GLfloat negRootT cyIRay dirVect x Val negRootPtr gt cylRay initPt yVal GLfloat negRootT cylRay dirVect yVal negRootPtr gt zVal cylRay initPt zVal GLfloat negRootT cylRay dirVect zVal return numRoots If poth points lie outside of the same cut plane then they do not represent a real intersection If only one lies outside then that ray needs to be intersected with the plane to find the true intersection with the cut cylinder If both points lie outside of different cut planes then both rays need to be intersected with their respective cut planes If neither lies outside the planes then they are the correct intersection points The cut planes are defined in scene coords The points need to be also int MovePointsWithinCutPlanes POINT3 posRootPtr POINT3 negRootPtr PLANE cutPl1 PLANE int goodPts lt TRUE GLdouble pt1PI2 21011 pt2PI2 RAY scrRay Substitute point value into planar eqn to determine which side of the plane the point is on 211111 cutPI1 xParam posRootPtr gt xVal cutPI 1 yParam posRootPtr gt yVal cutP11 zParam posRootPtr gt zVal
79. ter cutting the end put cylinder back so that the center of the bottom is at origin Calculate top cut plane At top extend top point up Z axis by enough that top cut plane will cut entirely through cylinder MoveToOrigin bot center AlignZAxis bot center top center Now the cylinder we want to build is aligned along the positive Z axis and its bottom is at the origin Move cylinder up above origin along Z axis enough to allow cut ExtendCyl angZ1 bot amp deltaZBot How much to move Now move origin down so cylinder will be longer by deltaZBot glTranslatef GLfloat 0 0 GLfloat 0 0 deltaZBot Extend top of cylinder ExtendCyl angZ2 top amp deltaZTop How much to extend newTop center zVal DistBetweenPoints top center bot center deltaZTop deltaZBot newTop center xVal top center x Val newTop center yVal top center y Val newTop radius top radius Create the cylinder extended so that the cuts can be made gluCylinder cylQuadric bot radius top radius newTop center zVal 10 10 Turn off clip planes so they don t affect any other objects glDisable GL CLIP PLANEO glDisable CLIP PLANEI density top density bot density 2 0 put ray data in pixel array MapCylinderPixels bot radius top radius newTop center zVal cutPl1 cutPI2 density bot temperature top temperature glPopMatrix Restore matrix
80. th one hundred images on a time scale measured in hours not days or weeks Because SXI is intended for use by solar physicists at perhaps several institutions it is written in a combination of the C language and OpenGL Both C and OpenGL are available on many platforms As a result SXI runs under both Windows95 and Unix SXI INTERNALS The actual code for SXI is included as Appendices A F Some of the more interesting aspects of the operation of SXI are discussed below Building the Plasma Loop The loop is constructed from common building blocks Since the radius of the loop can change along the length of the loop tapered cylindrical shapes are used Each cylinder is truncated at both ends to fit neatly against the next cylinder in the loop This process requires many steps The function MakeCylinder takes two DISKs each containing a center position and radius to mark the bottom and top of the cylinder A DISK is a structure defined in the file sxi h that is made up of the elements center of type POINT3 radius density and temperature A POINTS is a structure made up of three values for the x y and z coordinates of the point MakeCylinder also takes two more DISKs containing the prior and subsequent points and radii These extra DISKs are used to determine the planes which will cut the ends of the cylinder so that it will fit neatly to the next one in line It makes a tapered cylinder using the DISKs for the bottom
81. tic x ray image must be drawn in screen coordinates but the individual cylinders are initially divi symmetric about the z axis and their equations of form are relatively simple in that coordinate system hereafter called the cylinder coordinate system The rays emanating from the screen are known in screen coordinates but must be intersected with the cylinder in cylinder coordinates Thus the screen ray must be transformed to cylinder coordinates Since this transformation involves only translation and rotation any distances will remain unchanged Once the screen ray is in cylinder coordinates it must then be tested to see if it intersects with the truncated tapered cylinder If it does intersect then the distance through the cylinder is calculated The OpenGL Programming Guide Woo was the primary reference used for the OpenGL specific functions and for the various coordinate transformations The steps proceed as follows First calculate the direction vector from the screen into the scene Since we are using parallel projection there is only one direction vector for the entire scene Second the bounding rectangle for the cylinder in screen coordinates Foley 462 465 This is done by placing four points around the ends of the un clipped cylinder The positions of the four points are at square root of 2 times the radius of the end of the cylinder along both the x and y axes at z equals 0 and at z equals the cylinder height The sq
82. to solve the simultaneous equation for its roots which are the intersections We need the equation for the cylinder in question in cylinder coordinates this is the quadric The equation to be solved is in the form of k2 t k1 t 0 where t is the variable that determines how far along the ray we WC 80 to find an intersection For the case of a right tapered cylinder called cylinder for short k2 ax ax Q11 ay ay Q22 az az Q33 2 ax x0 Q11 2 ay y0 Q22 2 az z0 Q33 034 x0 x0 Q11 y0 y0 Q22 z0 z0 Q33 Q44 2 20 034 The Q s are non zero elements of the quadric form of a right tapered cylinder ax ay and az are the direction vector elements for the ray and x0 0 and 20 are the starting end point values for the ray 14 Image Storage 256 color lookup table is defined for the synthetic x ray image colors Each pixel color is stored as one 8 bit byte The value of the pixel byte is used to index the lookup table The pixel values are stored in an array with the lower left pixel in array position 0 This pixel array is exported as a bitmap file for use by other programs An example of a bitmap file converted to grayscale is shown in Figure 4 Figure 4 Grayscale Version of Synthetic X Ray Image 15 CONCLUSION AND POSSIBLE IMPROVEMENTS The SXI program meets its initial requirements It draws a synthetic solar plasma loop from the output of a mathe
83. tr gt xVal scrRay dirVect yVal negRootPtr gt yVal posRootPtr gt yVal scrRay dirVect zVal negRootPtr gt zVal posRootPtr gt zVal lt 0 if DEBUG_CLIP_PLANE cout lt lt lt lt endl DEBUG CLIP PLANE 0 IntersectRayWithCutPlane posRootPtr cutPl1 scrRay else if pt1P12 gt 0 if DEBUG CLIP PLANE cout lt lt pt1PI2 lt lt endl DEBUG CLIP PLANE 0 IntersectRayWithCutPlane posRootPtr cutP12 scrRay DEBUG CLIP PLANE 0 50 22211 lt 0 if DEBUG CLIP PLANE cout lt lt pt2PI1 lt lt endl DEBUG CLIP PLANE 0 IntersectRayWithCutPlane negRootPtr scrRay j else if pt2PI2 gt 0 if DEBUG CLIP PLANE lt lt pt2PI2 lt lt endl DEBUG CLIP PLANE 0 AC JM AUR negRootPtr cutPI2 scrRay return goodPts Using planar eqn Cz 4 0 and a ray defined by a direction vector and an initial point with the ray written in parametric form as x initPt xVal t dirVect xVal and similarly for y and 2 we can find the value of t at the point where the ray intersects the plane if there is an intersection From Foley VanDam Feiner Hughes and Phillips p 460 461 t A initPT xVal B initPt yVal C initPt zVal D A dirVect xVal B dirVect yVal C dirVect zVal If the denominator is 0 then the ray is parallel with the plane and the
84. ty sensitivity2 2 exposure 16 AdjustPixelColor xPos yPos counts Function to find the min and max of several bounding points void MinMax GLdouble minMax 20113 scrPt int init if init TRUE minMax 0 scrPt xVal minMax 1 scrPt xVal minMax 2 scrPt yVal minMax 3 scrPt y Val init FALSE else 44 if scrPt xVal lt minMax 0 minMax 0 scrPt xVal else if scrPt xVal gt minMax 1 minMax 1 7 scrPt xVal if scrPt yVal lt minMax 2 minMax 2 scrPt yVal else if scrPt yVal gt minMax 3 minMax 3 scrPt yVal He keke ke ok ke ok ok ok koe ke ok ke ke ke oe ke ke ke ke e ke he ok ke ok ke ok ok oe ok ke ok ke ook ke ok ke ok ke oe ok he oj ke Part 3 and Part 4 See if ray from this pixel intersects cylinder Make sure that ray intersects twice or else is not used Part 5 If ray intersects cylinder determine the distance between in and out points Equations for right truncated cone From z 0 to 2 20 2 0 have circle y y 1 1 At z 20 have circle x x y y 2 2 circle at zis x x y y r1 r1 r2 z zr2 2 or x x y y rl r1 12 z zr2 2 0 x x y y r1 r1 2r1 r1 r2 z zr2 r1 12 2 z zr2 2 0 x x y y r1 r1 2r1 r1 r2 z zr2 r1 12 2 z zr2 2 0 coeffx x 1 so 011 1 coeff y y 1 so
85. uare root of 2 times the radius makes certain that no matter what viewing angle about the z axis is used the points lie outside of the image of the cylinder end These eight points are transformed to screen pixel coordinates and the maximum and 11 minimum x and y values among all eight points are determined These maximum and minimum x and y values then determine a bounding box for the image of the cylinder in screen coordinates see Figure 2 As constructed Screen view with bounding box shown Figure 2 Creating Bounding Box for Cylinder Third for each pixel in the bounding rectangle determine if the ray emanating from it intersects with the cylinder This requires mapping the screen pixels and rays to the cylinder coordinates then calculating the intersections Fourth if the ray intersection with the cylinder is outside one of the clip planes for the ends of this cylinder replace this intersection with the intersection of the ray with the clip plane But if the second ray intersection with the cylinder is also outside the same clip plane then this ray doesn t intersect the clipped cylinder at all see Figure 3 12 Ray from screen pixel Intersection with back of cylinder Intersection with cut plane Original intersection Cylinder end plane Figure 3 Ray Intersection with Cut Plane Fifth calculate the distance between the two points of intersection Sixth set the x ray image bu

31762102746607 - ScholarWorks

Contents

Download Pdf Manuals

Related Search

Related Contents