123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342 |
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- #include <iostream>
- #include <fstream>
-
- #include "dipolesource.h"
- #include "layeredearth.h"
- #include "receiverpoints.h"
- #include "emearth1d.h"
- #include "WireAntenna.h"
- #include "PolygonalWireAntenna.h"
-
- #if LEMMAUSEVTK
- #include "vtkRenderer.h"
- #include "vtkRenderWindow.h"
- #include "vtkRenderWindowInteractor.h"
- #include "vtkRenderLargeImage.h"
- #include "vtkPNGWriter.h"
- #include "vtkRectilinearGrid.h"
- #include "vtkXMLRectilinearGridWriter.h"
- #include "vtkDoubleArray.h"
- #include "vtkFieldData.h"
- #include "vtkCellData.h"
- #endif
-
-
-
-
- #include "timer.h"
-
- using namespace Lemma;
-
- double randDouble(double low, double high) {
-
- double temp;
-
- if (low > high) {
- temp = low;
- low = high;
- high = temp;
- }
-
- temp = (rand() / (static_cast<double>(RAND_MAX) + 1.0))
- * (high - low) + low;
- return temp;
- }
-
-
- int main() {
-
-
- jsw_timer timer;
-
-
- srand(time(0));
-
- PolygonalWireAntenna *pa = PolygonalWireAntenna::New();
- pa->SetNumberOfFrequencies(1);
- pa->SetFrequency(0, 540000);
- Real Depth = 370;
- Real Width = 0.14;
- Real Height = 0.70;
-
-
- pa->SetNumberOfPoints(5);
-
-
-
- pa->SetPoint(0, Vector3r( -Width/2., -Height/2., Depth));
- pa->SetPoint(1, Vector3r( Width/2., -Height/2., Depth));
- pa->SetPoint(2, Vector3r( Width/2., Height/2., Depth));
- pa->SetPoint(3, Vector3r( -Width/2., Height/2., Depth));
- pa->SetPoint(4, Vector3r( -Width/2., -Height/2., Depth));
-
- pa->SetCurrent(1.);
- pa->SetNumberOfTurns(6);
-
-
-
-
-
-
-
-
-
-
-
- Real Sigma = 1.0/20. ;
- VectorXcr sigma(2);
- sigma << Complex(0.,0), Complex(Sigma,0);
- VectorXr thick(1);
- thick << 10;
- LayeredEarthEM *earth = LayeredEarthEM::New();
- earth->SetNumberOfLayers(2);
- earth->SetLayerConductivity(sigma);
-
-
-
- ReceiverPoints *receivers = ReceiverPoints::New();
- Vector3r loc;
- Real ox = -5.*Width - .003373;
- Real oy = -3.*Height - .003373;
- Real oz = Depth - Height/2. - .003373;
- int nx = 120;
- int ny = 180;
- int nz = 100;
- Real hx = 11.*Width/nx;
- Real hy = 7.*Height/ny;
- Real hz = 1.*Height/nz;
-
- receivers->SetNumberOfReceivers(nx*ny*nz);
- int ir = 0;
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
- loc << ox+ix*hx, oy+iy*hy, oz+iz*hz;
- receivers->SetLocation(ir, loc);
- ++ ir;
- }
- }
- }
-
-
- EMEarth1D *EmEarth = EMEarth1D::New();
-
- EmEarth->AttachWireAntenna(pa);
- EmEarth->AttachLayeredEarthEM(earth);
- EmEarth->AttachReceiverPoints(receivers);
- EmEarth->SetFieldsToCalculate(E);
-
- EmEarth->SetHankelTransformMethod(ANDERSON801);
-
-
- timer.begin();
- EmEarth->CalculateWireAntennaFields();
- Real paTime = timer.end();
- std::cout << "Polygonal wire antennae time: " << paTime << "\n";
-
-
-
-
-
-
-
-
-
- #if LEMMAUSEVTK
-
-
- vtkDoubleArray *xCoords = vtkDoubleArray::New();
- xCoords->InsertNextValue(ox-hx/2.);
- double xm1 = ox-hx/2.;
- for (int ix=0; ix<nx; ix++) {
- xCoords->InsertNextValue(xm1 + hx);
- xm1 += hx;
- }
-
- vtkDoubleArray *yCoords = vtkDoubleArray::New();
- yCoords->InsertNextValue(oy-hy/2.);
- double ym1 = oy-hy/2.;
- for (int iy=0; iy<ny; iy++) {
- yCoords->InsertNextValue(ym1 + hy);
- ym1 += hy;
- }
-
- vtkDoubleArray *zCoords = vtkDoubleArray::New();
- zCoords->InsertNextValue(oz-hz/2.);
- double zm1 = oz-hz/2.;
- for (int iz=0; iz<nz; iz++) {
- zCoords->InsertNextValue(zm1 + hz);
- zm1 += hz;
- }
-
- vtkDoubleArray *EReal = vtkDoubleArray::New();
- vtkDoubleArray *EImag = vtkDoubleArray::New();
- vtkDoubleArray *Watts = vtkDoubleArray::New();
- EReal->SetNumberOfComponents(3);
- EImag->SetNumberOfComponents(3);
- Watts->SetNumberOfComponents(1);
-
- ir = 0;
- Real WattsTotal(0);
- Real WattsInterior(0);
- for (int iz=0; iz<nz; ++iz) {
- for (int iy=0; iy<ny; ++iy) {
- for (int ix=0; ix<nx; ++ix) {
-
- Vector3cr E = receivers->GetEfield(0, ir);
- EReal-> InsertTuple3(ir, real(E(0)), real(E(1)), real(E(2)));
- EImag-> InsertTuple3(ir, imag(E(0)), imag(E(1)), imag(E(2)));
-
-
- if ( std::sqrt( std::pow(ox+ix*hx,2) + std::pow(oz+iz*hz - Depth, 2) ) < .085 && std::abs(oy+iy*hy) < Height/2. ) {
- Watts-> InsertTuple1(ir, 1e-20 );
- WattsInterior += .5* (( pow((std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
- } else {
- Watts-> InsertTuple1(ir, .5* (pow( (std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
-
- WattsTotal += .5* (( pow((std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
- }
- ++ ir;
- }
- }
- }
- std::cout << "Total Power: " << Sigma << "\t" << WattsTotal << "\t" << WattsInterior << endl;
-
- EReal->SetName("E_real");
- EImag->SetName("E_imag");
- Watts->SetName("Power");
-
- vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
- rgrid->SetDimensions(nx+1,ny+1,nz+1);
- rgrid->SetXCoordinates(xCoords);
- rgrid->SetYCoordinates(yCoords);
- rgrid->SetZCoordinates(zCoords);
-
- rgrid->GetCellData()->AddArray(EReal);
- rgrid->GetCellData()->AddArray(EImag);
- rgrid->GetCellData()->AddArray(Watts);
-
-
- vtkXMLRectilinearGridWriter *gridWrite = vtkXMLRectilinearGridWriter::New();
- gridWrite->SetInputData(rgrid);
- gridWrite->SetFileName("ors.vtr");
- gridWrite->Write();
-
-
- #endif
- EmEarth->Delete();
- receivers->Delete();
- earth->Delete();
-
- #if LEMMAUSEVTKX
-
- vtkRenderer *renderer = vtkRenderer::New();
- vtkRenderWindow *renWin = vtkRenderWindow::New();
- vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
-
- std::cout << "wire antennae approximating " << std::endl;
- loc << 0, .5*Width, Depth+.5*Height;
- pa->ApproximateWithElectricDipoles(loc);
- std::cout << "Wire approximated with " << pa->GetNumberOfDipoles() << std::endl;
-
- vtkActor **pdipActors = new vtkActor*[pa->GetNumberOfDipoles()];
- for (int id=0; id<pa->GetNumberOfDipoles(); ++id) {
- pdipActors[id] = pa->GetVtkActor(id);
- renderer->AddActor(pdipActors[id]);
- }
-
-
-
-
- renderer->SetBackground(1,1,1);
-
-
- renWin->AddRenderer(renderer);
- renWin->SetWindowName("Wire antennae");
-
- iren->SetRenderWindow(renWin);
- iren->Initialize();
- iren->Start();
- iren->Render();
-
- #if 0
- cout << "Enter File name?: ";
- std::string pngName;
- std::cin >> pngName;
- vtkPNGWriter *pngwrite = vtkPNGWriter::New();
- vtkRenderLargeImage *renlarge = vtkRenderLargeImage::New();
- renlarge->SetInput(renderer);
- renlarge->SetMagnification(2);
- pngwrite->SetInputConnection(renlarge->GetOutputPort());
- pngName.append(".png");
- pngwrite->SetFileName(pngName.c_str());
- pngwrite->Write();
- #endif
-
- #endif
-
-
-
-
-
- return 0;
- }
|