/source/pressureStar.c

https://github.com/EyNuel/cTraceo · C · 147 lines · 58 code · 19 blank · 70 comment · 7 complexity · 383b1eeb216ae404b81bd4cb4d784a2e MD5 · raw file

  1. /****************************************************************************************
  2. * pressureStar.c *
  3. * (formerly "star.for") *
  4. * Determines the star pressure contributions for particle velocity components of a *
  5. * specific ray at a specify coordinate. *
  6. * *
  7. * ------------------------------------------------------------------------------------ *
  8. * Website: *
  9. * https://github.com/EyNuel/cTraceo/wiki *
  10. * *
  11. * License: This file is part of the cTraceo Raytracing Model and is released under the *
  12. * Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License *
  13. * http://creativecommons.org/licenses/by-nc-sa/3.0/ *
  14. * *
  15. * NOTE: cTraceo is research code under active development. *
  16. * The code may contain bugs and updates are possible in the future. *
  17. * *
  18. * Written for project SENSOCEAN by: *
  19. * Emanuel Ey *
  20. * emanuel.ey@gmail.com *
  21. * Copyright (C) 2011 - 2013 *
  22. * Signal Processing Laboratory *
  23. * Universidade do Algarve *
  24. * *
  25. * cTraceo is the C port of the FORTRAN 77 TRACEO code written by: *
  26. * Orlando Camargo Rodriguez: *
  27. * Copyright (C) 2010 *
  28. * Orlando Camargo Rodriguez *
  29. * orodrig@ualg.pt *
  30. * Universidade do Algarve *
  31. * Physics Department *
  32. * Signal Processing Laboratory *
  33. * *
  34. * ------------------------------------------------------------------------------------ *
  35. * Inputs: *
  36. * settings: Pointer to structure containing all input info. *
  37. * ray: Pointer the the ray, who's influence we are determining. *
  38. * rHyd: Range coordinate of the hydrophone. *
  39. * zHyd: Depth coordinate of the hydrophone. *
  40. * q0: Value of q at beginning of ray. *
  41. * *
  42. * Outputs: *
  43. * pressure_H: A 3 element array containing the horizontal *
  44. * pressure components ( LEFT, CENTER, RIGHT). *
  45. * pressure_V: A 3 element array containing the vertical *
  46. * pressure components ( TOP, CENTER, BOTTOM). *
  47. * NOTE: there is redundacy in the above arrays: the center *
  48. * component is repeated for performance reasons. *
  49. * *
  50. * Return Value: *
  51. * 0: At least one of the star coordinates (rLeft, rRight, zBottom, *
  52. * zTop) is outside the rays' range or rBox. *
  53. * Do not use the outputs -doing so will result in a segfault *
  54. * 1: Outputs are valid. *
  55. * *
  56. ****************************************************************************************/
  57. #pragma once
  58. #include "globals.h"
  59. #include <complex.h>
  60. #include "bracket.c"
  61. uintptr_t pressureStar(settings_t*, ray_t*, double, double, double, complex double*, complex double[]);
  62. uintptr_t pressureStar( settings_t* settings, ray_t* ray, double rHyd, double zHyd, double q0, complex double* pressure_H, complex double* pressure_V){
  63. double rLeft, rRight, zTop, zBottom;
  64. uintptr_t iHyd;
  65. double dzdr, tauRay, zRay, qRay, width;
  66. complex double ampRay;
  67. /* start with determining the coordinates for which we will need to calculate
  68. * acoustic pressure.
  69. */
  70. rLeft = rHyd - settings->output.dr;
  71. rRight = rHyd + settings->output.dr;
  72. zBottom = zHyd - settings->output.dz;
  73. zTop = zHyd + settings->output.dz;
  74. DEBUG(9,"rL: %.2lf; rR: %.2lf; zD: %.2lf; zU: %.2lf;\n", rLeft, rRight, zBottom, zTop);
  75. //we have to check that these points are within the rBox:
  76. if( rLeft < settings->source.rbox1 ||
  77. rRight >= settings->source.rbox2 ){
  78. DEBUG(8, "Either rLeft or rRight are outside rBox\n");
  79. return 0;
  80. }
  81. //find out at what index of the ray coordinates the hydrophone is located:
  82. if( bracket(ray->nCoords, ray->r, rLeft, &iHyd) ){
  83. // NOTE: this block will not be run if the index returned by bracket() is out of bounds.
  84. DEBUG(8, "iHyd: %u => iRefl: %u\n", (uint32_t)iHyd, (uint32_t)ray->iRefl[iHyd]);
  85. if ( iHyd<ray->nCoords-1 ){
  86. getRayParameters(ray, iHyd, q0, rLeft, &dzdr, &tauRay, &zRay, &ampRay, &qRay, &width);
  87. DEBUG(8, "dzdr: %e, tauRay: %e, zRay: %e, ampRay: %e, qRay: %e, w: %e\n", dzdr, tauRay, zRay, cabs(ampRay), qRay, width);
  88. getRayPressureExplicit(settings, ray, iHyd, zHyd, tauRay, zRay, dzdr, ampRay, width, &pressure_H[LEFT]);
  89. }
  90. }else{
  91. DEBUG(6, "rLeft (%lf) can not be bracketed.\n", rLeft);
  92. return 0;
  93. }
  94. if( bracket(ray->nCoords, ray->r, rHyd, &iHyd)){
  95. /* NOTE: this block will not be run if the index returned by
  96. * bracket() is out of bounds.
  97. */
  98. DEBUG(8, "iHyd: %u => iRefl: %u\n", (uint32_t)iHyd, (uint32_t)ray->iRefl[iHyd]);
  99. if ( iHyd<ray->nCoords-1 ){
  100. getRayParameters(ray, iHyd, q0, rHyd, &dzdr, &tauRay, &zRay, &ampRay, &qRay, &width);
  101. DEBUG(8, "dzdr: %e, tauRay: %e, zRay: %e, ampRay: %e, qRay: %e, w: %e\n", dzdr, tauRay, zRay, cabs(ampRay), qRay, width);
  102. getRayPressureExplicit(settings, ray, iHyd, zTop, tauRay, zRay, dzdr, ampRay, width, &pressure_V[TOP]);
  103. getRayPressureExplicit(settings, ray, iHyd, zHyd, tauRay, zRay, dzdr, ampRay, width, &pressure_V[CENTER]);
  104. getRayPressureExplicit(settings, ray, iHyd, zBottom,tauRay, zRay, dzdr, ampRay, width, &pressure_V[BOTTOM]);
  105. pressure_H[CENTER] = pressure_V[CENTER];
  106. }
  107. }else{
  108. DEBUG(6, "rHyd (%lf) can not be bracketed.\n", rHyd);
  109. return 0;
  110. }
  111. if( bracket(ray->nCoords, ray->r, rRight, &iHyd)){
  112. /* NOTE: this block will not be run if the index returned by
  113. * bracket() is out of bounds.
  114. */
  115. DEBUG(8, "iHyd: %u => iRefl: %u\n", (uint32_t)iHyd, (uint32_t)ray->iRefl[iHyd]);
  116. if ( iHyd<ray->nCoords-1 ){
  117. //DEBUG(3, "r: %lf, z: %lf, amp:%lf, iHyd: %u\n", ray->r[iHyd], ray->z[iHyd], cabs(ray->amp[iHyd]), (uint32_t)iHyd);
  118. getRayParameters(ray, iHyd, q0, rRight, &dzdr, &tauRay, &zRay, &ampRay, &qRay, &width);
  119. DEBUG(8, "dzdr: %e, tauRay: %e, zRay: %e, ampRay: %e, qRay: %e, w: %e\n", dzdr, tauRay, zRay, cabs(ampRay), qRay, width);
  120. getRayPressureExplicit(settings, ray, iHyd, zHyd, tauRay, zRay, dzdr, ampRay, width, &pressure_H[RIGHT]);
  121. }
  122. }else{
  123. DEBUG(6, "rRight (%lf) can not be bracketed.\n", rRight);
  124. return 0;
  125. }
  126. //success
  127. return 1;
  128. }