<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>thSolver v1.0: /home/carlos/Work/Code+Data/Publicly Released Code/v1.0/dep/src/depolarization.c Source File</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
<link href="tabs.css" rel="stylesheet" type="text/css">
</head><body>
<!-- Generated by Doxygen 1.4.6 -->
<div class="tabs">
  <ul>
    <li><a href="index.html"><span>Main&nbsp;Page</span></a></li>
    <li id="current"><a href="files.html"><span>Files</span></a></li>
  </ul></div>
<div class="tabs">
  <ul>
    <li><a href="files.html"><span>File&nbsp;List</span></a></li>
    <li><a href="globals.html"><span>Globals</span></a></li>
  </ul></div>
<h1>/home/carlos/Work/Code+Data/Publicly Released Code/v1.0/dep/src/depolarization.c</h1><a href="depolarization_8c.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/******************************************************************************</span>
<a name="l00002"></a>00002 <span class="comment">* File      : depolarization.c                                                *</span>
<a name="l00003"></a>00003 <span class="comment">* Author    : Carlos Rosales Fernandez (carlos@ihpc.a-star.edu.sg)            *</span>
<a name="l00004"></a>00004 <span class="comment">* Date      : 01-09-2006                                                      *</span>
<a name="l00005"></a>00005 <span class="comment">* Revision  : 1.0                                                             *</span>
<a name="l00006"></a>00006 <span class="comment">*******************************************************************************</span>
<a name="l00007"></a>00007 <span class="comment">* DESCRIPTION                                                                 *</span>
<a name="l00008"></a>00008 <span class="comment">* Takes the semiaxis lengths of an ellipsoid in array axis, and calculates    *</span>
<a name="l00009"></a>00009 <span class="comment">* the depolarization factors numerically using a high quality Gauss-Legendre  *</span>
<a name="l00010"></a>00010 <span class="comment">* quadrature, returning them in array Ld.                                     *</span>
<a name="l00011"></a>00011 <span class="comment">******************************************************************************/</span>
<a name="l00012"></a>00012 
<a name="l00013"></a>00013 <span class="comment">/******************************************************************************</span>
<a name="l00014"></a>00014 <span class="comment">* COPYRIGHT &amp; LICENSE INFORMATION                                             *</span>
<a name="l00015"></a>00015 <span class="comment">*                                                                             *</span>
<a name="l00016"></a>00016 <span class="comment">* Copyright 2006 Carlos Rosales Fernandez and The Institute of High           *</span>
<a name="l00017"></a>00017 <span class="comment">* Performance Computing (A*STAR)                                              *</span>
<a name="l00018"></a>00018 <span class="comment">*                                                                             *</span>
<a name="l00019"></a>00019 <span class="comment">* This file is part of depSolver.                                             *</span>
<a name="l00020"></a>00020 <span class="comment">*                                                                             *</span>
<a name="l00021"></a>00021 <span class="comment">* depSolver is free software; you can redistribute it and/or modify           *</span>
<a name="l00022"></a>00022 <span class="comment">* it under the terms of the GNU General Public License as published by        *</span>
<a name="l00023"></a>00023 <span class="comment">* the Free Software Foundation; either version 2 of the License, or           *</span>
<a name="l00024"></a>00024 <span class="comment">* (at your option) any later version.                                         *</span>
<a name="l00025"></a>00025 <span class="comment">*                                                                             *</span>
<a name="l00026"></a>00026 <span class="comment">* depSolver is distributed in the hope that it will be useful,                *</span>
<a name="l00027"></a>00027 <span class="comment">* but WITHOUT ANY WARRANTY; without even the implied warranty of              *</span>
<a name="l00028"></a>00028 <span class="comment">* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               *</span>
<a name="l00029"></a>00029 <span class="comment">* GNU General Public License for more details.                                *</span>
<a name="l00030"></a>00030 <span class="comment">*                                                                             *</span>
<a name="l00031"></a>00031 <span class="comment">* You should have received a copy of the GNU General Public License           *</span>
<a name="l00032"></a>00032 <span class="comment">* along with depSolver; if not, write to the Free Software                    *</span>
<a name="l00033"></a>00033 <span class="comment">* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA  *</span>
<a name="l00034"></a>00034 <span class="comment">******************************************************************************/</span>
<a name="l00035"></a>00035 
<a name="l00036"></a>00036 <span class="preprocessor">#include "<a class="code" href="depolarization_8h.html">depolarization.h</a>"</span>
<a name="l00037"></a>00037 
<a name="l00038"></a><a class="code" href="forceEllipsoid__tria6_8h.html#9585e7af1008faa928783689ecaab150">00038</a> <span class="keywordtype">int</span> <a class="code" href="depolarization_8c.html#9585e7af1008faa928783689ecaab150">depolarization</a>(<span class="keywordtype">double</span> *axis, <span class="keywordtype">double</span> *Ld)
<a name="l00039"></a>00039 {
<a name="l00040"></a>00040     <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span>  i;
<a name="l00041"></a>00041     <span class="keywordtype">double</span>  A, a2, b2, c2, J, L, R, W, x;
<a name="l00042"></a>00042     <span class="keywordtype">double</span>  N[2], XL[2];
<a name="l00043"></a>00043     <span class="keywordtype">double</span>  X[2][2];
<a name="l00044"></a>00044 
<a name="l00045"></a>00045     <span class="comment">/* Initialize */</span>
<a name="l00046"></a>00046     Ld[0] = Ld[1] = Ld[2] = 0.0;
<a name="l00047"></a>00047     X[0][0] = X[0][1] = X[1][1] = 0.0;
<a name="l00048"></a>00048     A = 0.5*axis[0]*axis[1]*axis[2];
<a name="l00049"></a>00049     a2 = axis[0]*axis[0];
<a name="l00050"></a>00050     b2 = axis[1]*axis[1];
<a name="l00051"></a>00051     c2 = axis[2]*axis[2];
<a name="l00052"></a>00052 
<a name="l00053"></a>00053     <span class="comment">/* Integrate for Lx using Gaussian quadrature */</span>
<a name="l00054"></a>00054     X[1][0] = 100*a2;
<a name="l00055"></a>00055     <span class="keywordflow">for</span>(i = 0; i &lt; 64; i ++){
<a name="l00056"></a>00056         L = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][0];
<a name="l00057"></a>00057         W = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][1];
<a name="l00058"></a>00058     
<a name="l00059"></a>00059         <a class="code" href="depolarization_8h.html#ec879a354e79715bcf18d7043627779e">shape_line2</a>(L,N);
<a name="l00060"></a>00060         J = <a class="code" href="depolarization_8h.html#435064417821bf36371611d6c34ef0d7">X2L_line2</a>(X,XL,N);
<a name="l00061"></a>00061         x = XL[0];
<a name="l00062"></a>00062         R = J*W/sqrt((x+a2)*(x+b2)*(x+c2));
<a name="l00063"></a>00063     
<a name="l00064"></a>00064         Ld[0] += R/(x+a2);
<a name="l00065"></a>00065     }
<a name="l00066"></a>00066 
<a name="l00067"></a>00067     <span class="comment">/* Integrate for Ly using Gaussian quadrature */</span>
<a name="l00068"></a>00068     X[1][0] = 100*b2;
<a name="l00069"></a>00069     <span class="keywordflow">for</span>(i = 0; i &lt; 64; i ++){
<a name="l00070"></a>00070         L = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][0];
<a name="l00071"></a>00071         W = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][1];
<a name="l00072"></a>00072     
<a name="l00073"></a>00073         <a class="code" href="depolarization_8h.html#ec879a354e79715bcf18d7043627779e">shape_line2</a>(L,N);
<a name="l00074"></a>00074         J = <a class="code" href="depolarization_8h.html#435064417821bf36371611d6c34ef0d7">X2L_line2</a>(X,XL,N);
<a name="l00075"></a>00075         x = XL[0];
<a name="l00076"></a>00076         R = J*W/sqrt((x+a2)*(x+b2)*(x+c2));
<a name="l00077"></a>00077     
<a name="l00078"></a>00078         Ld[1] += R/(x+b2);
<a name="l00079"></a>00079     }
<a name="l00080"></a>00080 
<a name="l00081"></a>00081     <span class="comment">/* Integrate for Lz using Gaussian quadrature */</span>
<a name="l00082"></a>00082     X[1][0] = 100*c2;
<a name="l00083"></a>00083     <span class="keywordflow">for</span>(i = 0; i &lt; 64; i ++){
<a name="l00084"></a>00084         L = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][0];
<a name="l00085"></a>00085         W = <a class="code" href="depolarization_8h.html#779f2ad293a6b0f3d72fa6b0346d3e6a">GaussH</a>[i][1];
<a name="l00086"></a>00086     
<a name="l00087"></a>00087         <a class="code" href="depolarization_8h.html#ec879a354e79715bcf18d7043627779e">shape_line2</a>(L,N);
<a name="l00088"></a>00088         J = <a class="code" href="depolarization_8h.html#435064417821bf36371611d6c34ef0d7">X2L_line2</a>(X,XL,N);
<a name="l00089"></a>00089         x = XL[0];
<a name="l00090"></a>00090         R = J*W/sqrt((x+a2)*(x+b2)*(x+c2));
<a name="l00091"></a>00091     
<a name="l00092"></a>00092         Ld[2] += R/(x+c2);
<a name="l00093"></a>00093     }
<a name="l00094"></a>00094 
<a name="l00095"></a>00095 
<a name="l00096"></a>00096     Ld[0] = A*Ld[0];
<a name="l00097"></a>00097     Ld[1] = A*Ld[1];
<a name="l00098"></a>00098     Ld[2] = A*Ld[2];
<a name="l00099"></a>00099 
<a name="l00100"></a>00100     <span class="comment">/* Safety checks */</span>
<a name="l00101"></a>00101     <span class="keywordflow">for</span>(i = 0; i &lt; 3; i++){
<a name="l00102"></a>00102         <span class="keywordflow">if</span>(Ld[i] &lt; 0.0) <a class="code" href="depolarization_8h.html#2e4b1df9cfea8f7b71137f266c6df7a1">errorHandler</a>(<span class="stringliteral">"Error: depolarization() - factors must be positive!"</span>);
<a name="l00103"></a>00103     }
<a name="l00104"></a>00104     <span class="keywordflow">if</span>(fabs(1.0 - Ld[0] - Ld[1] - Ld[2]) &gt; 0.05)
<a name="l00105"></a>00105         printf(<span class="stringliteral">"Warning: depolarization() - error is larger than 0.05"</span>);
<a name="l00106"></a>00106 
<a name="l00107"></a>00107     <span class="keywordflow">return</span> 0;
<a name="l00108"></a>00108 }
</pre></div><hr size="1"><address style="align: right;"><small>Generated on Fri Sep 1 18:53:06 2006 for thSolver v1.0 by&nbsp;
<a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.4.6 </small></address>
</body>
</html>
