From 9a4a8954639e7a206263d5d7c2413248dbed0249 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 3 Jun 2020 14:20:02 -0300 Subject: [PATCH 01/30] Modeling step with madagascarRecipes package #2 Run 'scons -f SConscript' to generate the model and datacubes using defined function of madagascarRecipes package --- experiments/README.md | 1 + experiments/velocityInversion/README.md | 1 + experiments/velocityInversion/SConscript | 64 ++ .../madagascarRecipes/LICENSE | 674 ++++++++++++++++++ .../madagascarRecipes/README.md | 2 + .../madagascarRecipes/__init__.py | 0 .../madagascarRecipes/errorReport.py | 54 ++ .../madagascarRecipes/generateSelfdoc | 27 + .../madagascarRecipes/kimodel.py | 187 +++++ .../madagascarRecipes/pefInterpolation.py | 156 ++++ .../madagascarRecipes/plot.py | 52 ++ .../madagascarRecipes/velocityAnalisys.py | 100 +++ 12 files changed, 1318 insertions(+) create mode 100644 experiments/README.md create mode 100644 experiments/velocityInversion/README.md create mode 100644 experiments/velocityInversion/SConscript create mode 100644 experiments/velocityInversion/madagascarRecipes/LICENSE create mode 100644 experiments/velocityInversion/madagascarRecipes/README.md create mode 100644 experiments/velocityInversion/madagascarRecipes/__init__.py create mode 100644 experiments/velocityInversion/madagascarRecipes/errorReport.py create mode 100644 experiments/velocityInversion/madagascarRecipes/generateSelfdoc create mode 100644 experiments/velocityInversion/madagascarRecipes/kimodel.py create mode 100644 experiments/velocityInversion/madagascarRecipes/pefInterpolation.py create mode 100644 experiments/velocityInversion/madagascarRecipes/plot.py create mode 100644 experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py diff --git a/experiments/README.md b/experiments/README.md new file mode 100644 index 0000000..5f21e89 --- /dev/null +++ b/experiments/README.md @@ -0,0 +1 @@ +#### Experiments and studies of CRE Velocity Inversion diff --git a/experiments/velocityInversion/README.md b/experiments/velocityInversion/README.md new file mode 100644 index 0000000..8b76cfb --- /dev/null +++ b/experiments/velocityInversion/README.md @@ -0,0 +1 @@ +#### Velocity Inversion algorithm diff --git a/experiments/velocityInversion/SConscript b/experiments/velocityInversion/SConscript new file mode 100644 index 0000000..819fa2f --- /dev/null +++ b/experiments/velocityInversion/SConscript @@ -0,0 +1,64 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# SConscript (Madagascar Script) +# +# Purpose: Generate model and data cubes. +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 03/06/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Madagascar package +from rsf.proj import * + +# Error report function +import atexit +from madagascarRecipes.errorReport import print_build_failures + +# CRE recipe +from madagascarRecipes.pefInterpolation import pefInterpolation as pefin +from madagascarRecipes.kimodel import multiLayerModel as mlmod +from madagascarRecipes.kimodel import kirchoffNewtonModeling as kinewmod + +xmax = 6.0 +zmax = 2.0 + +layers = ((0.30,0.50,0.20,0.30), + (1.65,1.85,1.55,1.65)) + +velocities = (1.508, + 1.690, + 2.0) + +# Generate multi layer model and data cube +mlmod(interfaces='interfaces', + dipsfile='interfacesDip', + modelfile='mod1', + xmax=xmax, + zmax=zmax, + layers=layers, + velocities=velocities) + +kinewmod(reflectors='interfaces', + reflectorsDip='interfacesDip', + filename='multiLayerDataCube', + velocities=velocities) + +# Use aliases to split building +Alias('model','multiLayerDataCube.rsf') + +# Show error message if fail + +message = ''' +SCosncript Failed build +''' + +atexit.register(print_build_failures,message) +End() diff --git a/experiments/velocityInversion/madagascarRecipes/LICENSE b/experiments/velocityInversion/madagascarRecipes/LICENSE new file mode 100644 index 0000000..f288702 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/experiments/velocityInversion/madagascarRecipes/README.md b/experiments/velocityInversion/madagascarRecipes/README.md new file mode 100644 index 0000000..9448f2b --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/README.md @@ -0,0 +1,2 @@ +# madagascarRecipes +Personal library of Madagascar recipes diff --git a/experiments/velocityInversion/madagascarRecipes/__init__.py b/experiments/velocityInversion/madagascarRecipes/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/experiments/velocityInversion/madagascarRecipes/errorReport.py b/experiments/velocityInversion/madagascarRecipes/errorReport.py new file mode 100644 index 0000000..8509563 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/errorReport.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# errorReport.py (Madagascar Recipe) +# +# Purpose: SCons error report function. +# +# Important!: It should be called from a SConstruct +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 01/06/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Selfdoc string +''' +Python function to report errors in SCons + +If build fail report specific message. Import this module and add the following lineis to your SConstruct: + +import atexit +from madagascarRecipes.errorReport import print_build_failures + +Call error message by the following: + +atexit.register(print_build_failures,message) + +It will call function print_build_failures to report message variable that should be a string with the error message to be reported +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +__author__="Rodolfo Dirack " +__version__="1.0" + +def print_build_failures(message): + ''' + Print error message when build fail + :param message: string, error message + ''' + from SCons.Script import GetBuildFailures + for bf in GetBuildFailures(): + print("BUILD ERROR REPORT".center(40,'*')) + print("FAILED BUILD: %s" % (bf.node)) + print("MESSAGE REPORT".center(40,'*')) + print("%s" % (message)) + diff --git a/experiments/velocityInversion/madagascarRecipes/generateSelfdoc b/experiments/velocityInversion/madagascarRecipes/generateSelfdoc new file mode 100644 index 0000000..4f98839 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/generateSelfdoc @@ -0,0 +1,27 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# generateSelfdoc (SConscript) +# +# Purpose: Generate selfdoc for Madagascar recipes using pydoc. +# +# Important!: It should be called from 'scons' with the command +# 'scons -f generateSelfdoc' +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 04/03/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +import pydoc + +RECIPES=['kimodel', + 'pefInterpolation'] + +for i in RECIPES: + pydoc.writedoc(i) diff --git a/experiments/velocityInversion/madagascarRecipes/kimodel.py b/experiments/velocityInversion/madagascarRecipes/kimodel.py new file mode 100644 index 0000000..177b625 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/kimodel.py @@ -0,0 +1,187 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# kimodel.py (Madagascar Recipe) +# +# Purpose: Recipe to Kirchhoff modeling. +# +# Important!: It should be called from a SConstruct +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 04/03/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Selfdoc string +''' +Madagascar recipe to Kirchhoff modeling + +Define functions to generate subsurface models, apply Kirchhoff modeling +and Kirchhoff Newton modeling. +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +# Madagascar package +from rsf.proj import * + +__author__="Rodolfo Dirack " +__version__="1.0" + +def arr2str(array,sep=' '): + ''' + Convert a tuple into a comma separeted string + ''' + return string.join(map(str,array),sep) + + +def multiLayerModel( + interfaces, + dipsfile, + modelfile, + xmax, + zmax, + layers, + velocities + ): + ''' + + Generate a multi layer model to use in the program sfkirmod_newton + + :out interfaces: RSF filename, interpolated interfaces file + :out dipsfile: RSF filename, dips of interfaces + :out modelfile: RSF filename, model for ploting + :param xmax: interger, max x axis model distance + :param zmax: interger, max z axis model depth + :param layers: tuple array, points describing interfaces + :param velocities: tuple, layers velocities in Km/s + ''' + vstr=arr2str(velocities,',') + + n1 = len(layers[0]) + n2 = len(layers) + + Flow('layers.asc',None, + ''' + echo %s + n1=%d n2=%d o1=0 d1=%g + data_format=ascii_float in=$TARGET + ''' % (string.join(map(arr2str,layers),' '), + n1,n2,xmax/(n1-1))) + Flow('layers','layers.asc','dd form=native') + + d = 0.0101 # non-round for reproducibility + + Flow(interfaces,'layers', + 'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d))) + Flow(dipsfile,interfaces,'deriv scale=y') + + Flow(modelfile,interfaces, + ''' + unif2 d1=%g n1=%d v00=%s + ''' % (d,int(1.5+zmax/d),vstr)) + +def kirchoffNewtonModeling( + reflectors, + reflectorsDip, + filename, + velocities): + ''' + Kirchhoff modeling for a multi layer model + ''' + vstr=arr2str(velocities,',') + + # Kirchoff modeling for multi layer model + Flow(filename,[reflectors,reflectorsDip], + ''' + kirmod_newton nt=1001 dt=0.004 freq=10 + ns=201 ds=0.025 nh=161 dh=0.025 h0=0 s0=0 verb=y cmp=y + vstatus=0 velocity=%s debug=n fwdxini=y + xref=0 zref=0 dip=${SOURCES[1]} | + put d2=0.0125 label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s + ''' % (vstr)) + +def gaussianReflector(filename='gaussianReflector'): + ''' + + Generate a gaussian reflector to use in the program sfkirmod + + :param filename: RSF filename, gaussian reflector file + ''' + + # Modeling: Gaussian reflector in a velocity linear model + # velocity increases with depth with a 0.5 velocity gradient + Flow('gaussianReflector',None, + ''' + math d1=0.01 n1=2001 o1=-5 unit1=km label1=Offset + output="4-3*exp(-(x1-5)^2/9)" + ''') + + # Velocity Model + Flow('velocityModel','gaussianReflector', + ''' + window min1=0 max1=10 | + spray axis=1 n=451 d=0.01 o=0 label=Depth unit=km | + math output="1.5+0.5*x1+0.0*x2" + ''') + + Flow('reflectorDip','gaussianReflector','math output="2/3*(x1-5)*input" ') + + +def kirchoffModeling( + filename, + reflector, + reflectorDip, + nh, + dh, + h0, + ns, + ds, + s0, + freq, + dt, + nt, + v0, + gradz + ): + ''' + + Kirchhoff modeling function. Generate (time,cmp,offset) datacube for + a given v(z) model (Velocity increases with depth). + + :out filename: RSF filename, Modeled data cube + :param reflector: RSF filename, interfaces file + :param reflectorDip: RSF filename, interfaces dips file + :param nh: integer, number of offsets + :param dh: float, offset sampling + :param h0: float, first offset + :param ns: integer, number of shots/cmps + :param ds: float, shots/cmps sampling + :param s0: float, first shot/cmp position + :param freq: integer, ricker pulse frequency + :param dt: float, time sampling + :param nt: float, number of time samples + :param v0: float, near surface velocity + :param gradz: float, velocity gradient + ''' + + halfOffset = dh/2. + + # Kirchoff Modeling + Flow(filename,[reflector, reflectorDip], + ''' + kirmod cmp=y dip=${SOURCES[1]} + nh=%i dh=%g h0=%g + ns=%i ds=%g s0=%g + freq=%i dt=%g nt=%i + vel=%g gradz=%g gradx=0.0 verb=y | + put d2=%g label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s + ''' % (nh,dh,h0,ns,ds,s0,freq,dt,nt,v0,gradz,halfOffset)) + diff --git a/experiments/velocityInversion/madagascarRecipes/pefInterpolation.py b/experiments/velocityInversion/madagascarRecipes/pefInterpolation.py new file mode 100644 index 0000000..495b635 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/pefInterpolation.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# pefInterpolation.py (Madagascar Recipe) +# +# Purpose: Recipe to Preditive Adaptative Error Filters (PEF) interpolation +# to double CMP number of samples. +# +# Important!: It should be called from a SConstruct +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 04/03/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Selfdoc string +''' +Madagascar recipe to PEF interpolation + +Interpolate and increase CMP sampling using Preditive Adaptative Error Filters (PEF) interpolation. +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +__author__="Rodolfo Dirack " +__version__="1.0" + +# Madagascar package +from rsf.proj import * + +def pefInterpolation( + dataCube, + interpolated, + nm, + dm, + nt, + dt, + nhi, + a1, + a2, + rect1, + rect2 + ): + ''' + + PEF interpolation to double CMP sampling of the data cube + + :param dataCube: RSF filename, Seismic data cube to interpolate + :param interpolated: RSF filename, Interpolated seismic data cube + :param nm: integer, number of CMPs in the seismic data cube + :param dm: float, CMP sampling + :param nt: integer, number of time samples + :param dt: float, time sampling + :param nhi: integer, number of constant offsets gathers to interpolate + :param a1: integer, Number of PEF coeficients in time axis + :param a2: integer, Number of PEF coeficients in space axis + :param rect1: integer, Smooth radius in time + :param rect2: integer, Smooth radius in space + ''' + + # Divide CMP sampling + dm = dm/2 + + # Define mask file names using input filename + mask1 = dataCube+'-mask1' + mask = dataCube+'-mask' + aa = dataCube+'-aa' + bb = dataCube+'-bb' + a = dataCube+'-a' + b = dataCube+'-b' + zeroTraceGather = dataCube+'-zeroedGather' + mask0 = dataCube+'-mask0' + + + # Build a mask to interleave zero traces with original data traces + Flow(aa,None,'spike n1=%i d1=%g o1=0' %(nm,dm)) + Flow(bb,None,'spike n1=%i d1=%g o1=0 mag=0' % (nm,dm)) + Flow(mask1,[bb, aa], + ''' + interleave axis=1 ${SOURCES[1]} | + dd type=int + ''') + + Flow(a,None,'spike n1=%i d1=%g o1=0' % (nm,dm)) + Flow(b,None,'spike n1=%i d1=%g o1=0 mag=0' % (nm,dm)) + Flow(mask,[a, b], + ''' + interleave axis=1 ${SOURCES[1]} | + dd type=int + ''') + Flow(zeroTraceGather,b, + ''' + spray axis=2 n=%i d=%g | + transp | + put label2=Offset unit2=Km label1=Time unit1=s + ''' %(nt,dm)) + + # Data Mask with double of traces in CMP (half of CMP sampling) + # Keep the same Time and Offset original data sampling + Flow(mask0,mask, + ''' + spray axis=1 n=%i d=%g + ''' %(nt,dt)) + + totalPefIterations = 100 + totalInterpolationIterations = 50 + + offsetGathers = [] + for offsetGatherIndex in range(nhi): + + offsetGather = dataCube+"-offsetGather-%i" % offsetGatherIndex + resampledOffsetGather = dataCube+"-resampledGather-%i" % offsetGatherIndex + interpolatedOffsetGather = dataCube+"-interpolatedGather-%i" % offsetGatherIndex + pefCoeficients = dataCube+"-pefCoeficients-%i" % offsetGatherIndex + + Flow(offsetGather,dataCube, + ''' + window n2=1 f2=%i + ''' % (offsetGatherIndex)) + + Flow(resampledOffsetGather,[offsetGather,zeroTraceGather], + ''' + interleave axis=2 ${SOURCES[1]} + ''') + + # Calculate adaptive PEF coeficients + Flow(pefCoeficients,[resampledOffsetGather,mask0], + ''' + apef jump=2 a=%i,%i rect1=%i rect2=%i niter=%g verb=y + maskin=${SOURCES[1]} + ''' % (a1,a2,rect1,rect2,totalPefIterations)) + + # Interpolation + Flow(interpolatedOffsetGather, [resampledOffsetGather,pefCoeficients,mask0,mask1], + ''' + miss4 exact=y filt=${SOURCES[1]} mask=${SOURCES[2]} niter=%g verb=y | + put d2=%g + ''' % (totalInterpolationIterations,dm)) + + offsetGathers.append(interpolatedOffsetGather) + + # Concatenate interpolated sections + Flow(interpolated,offsetGathers, + ''' + rcat axis=3 ${SOURCES[1:%d]} | + transp plane=23 + ''' % nhi) + + diff --git a/experiments/velocityInversion/madagascarRecipes/plot.py b/experiments/velocityInversion/madagascarRecipes/plot.py new file mode 100644 index 0000000..3715a31 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/plot.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# plot.py (Madagascar Recipe) +# +# Purpose: Recipe to generate plot files (VPL). +# +# Important!: It should be called from a SConstruct +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 04/03/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Selfdoc string +''' +Madagascar recipe to plot + +Define functions to generate vpl files. +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +# Madagascar package +from rsf.proj import * + +__author__="Rodolfo Dirack " +__version__="1.0" + +def velplot(title,label1='Depth',unit1='Km'): + ''' + Velocity ploting function + :param title: RSF filename, velocity file + :param label1: string, first axis label + :param unit1: string, first axis unit + ''' + + return ''' + grey color=j allpos=y title="%s" scalebar=y + barlabel=Velocity barunit=Km/s + label1="%s" unit1="%s" label2=Lateral unit2=Km + barreverse=y pclip=100 bias=2.8 + ''' % (title,label1,unit1) + + diff --git a/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py b/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py new file mode 100644 index 0000000..96a5826 --- /dev/null +++ b/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py @@ -0,0 +1,100 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# velocityAnalisys.py () +# +# Purpose: Madagascar recipe to do velocity analisys. +# +# Site: https://dirack.github.io +# +# Version 1.0 +# +# Programer: Rodolfo A C Neves (Dirack) 01/04/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +from rsf.proj import * +import math + +# Selfdoc string +''' +Madagascar recipe to Velocity analisys with automatic velocity scan and picking + +Define functions to do a velocity analisys of a (time,offset,CMP) datacube. +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +# Madagascar package +from rsf.proj import * + +__author__="Rodolfo Dirack " +__version__="1.0" + +def velocityAnalisys( + dataCube, + pick, + nmo, + vrms, + v0, + dv, + nv, + vel0, + rect1=15, + rect2=40, + rect3=3, + dt=0.004) + ''' + Velocity analisys and automatic picking and semblance + :param datacube: RSF filename, (time, offset, CMP) datacube + :out pick: RSF filename, picking velocity + :out stack: RSF filename, stacked section + :out vrms: RSF filename, vrms velocity + :param dt: float, time sampling + :param v0: float, inicial scan velocity + :param dv: float, scan velocity sampling + :param nv: int, number of scan velocities + :param vel0: float, initial velocity picking + :param rect1: int, smooth time radius to velocity scan + :param rect2: int, smooth space radius to velocity scan + :param rect3: int, smooth time radius post picking + ''' + + scan = 'scan'+dataCube + nmo = 'nmo'+dataCube + + # Velocity scan + Flow(scan,[dataCube], + ''' + vscan semblance=y + v0=%g dv=%g nv=%d + ''' % (v0,dv,nv)) + + # Velocity picking + Flow(pick,scan, + ''' + pick rect1=%i rect2=%i vel0=%g smooth=y | + smooth rect1=%i + ''' %(rect1,rect2,vel0,rect3)) + + # NMO and stack + Flow(nmo,[dataCube,pick], + 'nmo half=n velocity=${SOURCES[1]}') + + Flow(stack,nmo,'stack') + + # generate RMS velocity : + # sfmul multiply input data by itself to make v**2 + # sfcausint integrate v**2 (this does not include scaling by dt) + # sfmath scale by dt and divide by t + Flow(vrms,pick, + ''' + sfmul $SOURCE | + sfcausint | + sfmath output="sqrt(input*%g/(x1+%g))" | + sfput n3=1 d3=1 o3=0 + '''%(dt,dt)) From cf93671b3be65cfa205236bb78f9a8947eba3b18 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 3 Jun 2020 15:20:16 -0300 Subject: [PATCH 02/30] Ignore madagascar files --- experiments/velocityInversion/.gitignore | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 experiments/velocityInversion/.gitignore diff --git a/experiments/velocityInversion/.gitignore b/experiments/velocityInversion/.gitignore new file mode 100644 index 0000000..b0a5fac --- /dev/null +++ b/experiments/velocityInversion/.gitignore @@ -0,0 +1,4 @@ +*.rsf +*.vpl +*.asc +*.pyc From 14455c9af733d39ed97ebb7b2f8e102bf00c6758 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 3 Jun 2020 15:21:36 -0300 Subject: [PATCH 03/30] Correct bugs in velocityAnalisys recipe #2 This bug correction waas already reported to oficial madagascarRecipes repository and will be corrected in next version --- .../madagascarRecipes/velocityAnalisys.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py b/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py index 96a5826..3a03fed 100644 --- a/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py +++ b/experiments/velocityInversion/madagascarRecipes/velocityAnalisys.py @@ -35,10 +35,9 @@ __author__="Rodolfo Dirack " __version__="1.0" -def velocityAnalisys( - dataCube, +def velocityAnalisys(dataCube, pick, - nmo, + stack, vrms, v0, dv, @@ -47,7 +46,7 @@ def velocityAnalisys( rect1=15, rect2=40, rect3=3, - dt=0.004) + dt=0.004): ''' Velocity analisys and automatic picking and semblance :param datacube: RSF filename, (time, offset, CMP) datacube From 0d5388ca7acd654e814b2e794006d08a32390d8f Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 3 Jun 2020 15:22:57 -0300 Subject: [PATCH 04/30] Kirchhoff time migration #2 Generate Kirchhoff time migrated section after stacking --- experiments/velocityInversion/SConscript | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/experiments/velocityInversion/SConscript b/experiments/velocityInversion/SConscript index 819fa2f..82dc1f5 100644 --- a/experiments/velocityInversion/SConscript +++ b/experiments/velocityInversion/SConscript @@ -26,6 +26,7 @@ from madagascarRecipes.errorReport import print_build_failures from madagascarRecipes.pefInterpolation import pefInterpolation as pefin from madagascarRecipes.kimodel import multiLayerModel as mlmod from madagascarRecipes.kimodel import kirchoffNewtonModeling as kinewmod +from madagascarRecipes.velocityAnalisys import velocityAnalisys as nmoStack xmax = 6.0 zmax = 2.0 @@ -51,6 +52,25 @@ kinewmod(reflectors='interfaces', filename='multiLayerDataCube', velocities=velocities) +nmoStack(dataCube='multiLayerDataCube', + pick='vnmo', + stack='stackedSection', + vrms='vrms', + v0=1.5, + dv=0.01, + nv=100, + vel0=1.5, + rect1=15, + rect2=40, + rect3=3, + dt=0.004) + +Flow('kirchhoffSection',['stackedSection','vrms'], + ''' + kirchnew velocity=${SOURCES[1]} + ''' + ) + # Use aliases to split building Alias('model','multiLayerDataCube.rsf') From 4a1793ff775bd93d7ca74fe8e2e9df5654ef2c8d Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 4 Jun 2020 22:36:56 -0300 Subject: [PATCH 05/30] Time to depth conversion #2 Simple time to depth conversion after Kirchhoff Migration --- experiments/velocityInversion/SConscript | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/experiments/velocityInversion/SConscript b/experiments/velocityInversion/SConscript index 82dc1f5..6e261e1 100644 --- a/experiments/velocityInversion/SConscript +++ b/experiments/velocityInversion/SConscript @@ -3,7 +3,8 @@ # # SConscript (Madagascar Script) # -# Purpose: Generate model and data cubes. +# Purpose: Generate model and data cubes, NMO stacking, +# Kirchhoff Migration and time to depth conversion. # # Site: http://www.dirackslounge.online # @@ -52,6 +53,7 @@ kinewmod(reflectors='interfaces', filename='multiLayerDataCube', velocities=velocities) +# NMO Stack nmoStack(dataCube='multiLayerDataCube', pick='vnmo', stack='stackedSection', @@ -65,12 +67,22 @@ nmoStack(dataCube='multiLayerDataCube', rect3=3, dt=0.004) +# Kirchhoff Migration Flow('kirchhoffSection',['stackedSection','vrms'], ''' kirchnew velocity=${SOURCES[1]} ''' ) +# Time to depth Conversion +nz = 199 +dz = 0.01 +Flow('depthVelModel',['kirchhoffSection','vrms'], + ''' + time2depth velocity=${SOURCES[1]} nz=%i dz=%g intime=y | + put label1=Depth unit1=Km + ''' % (nz,dz)) + # Use aliases to split building Alias('model','multiLayerDataCube.rsf') From aa08585e8598bda7b1cb307828271bc4fb180f58 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Fri, 5 Jun 2020 17:15:55 -0300 Subject: [PATCH 06/30] Generate Depth velocity model and Depth section #2 Those sections are a by product of time to depth conversion --- experiments/velocityInversion/SConscript | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/experiments/velocityInversion/SConscript b/experiments/velocityInversion/SConscript index 6e261e1..992d4ee 100644 --- a/experiments/velocityInversion/SConscript +++ b/experiments/velocityInversion/SConscript @@ -74,20 +74,37 @@ Flow('kirchhoffSection',['stackedSection','vrms'], ''' ) +# Convert RMS velocity to Dix velocity +Flow(['vdix','vout'],'vrms','dix vrmsout=${TARGETS[1]}') + # Time to depth Conversion +# Depth velocity model nz = 199 dz = 0.01 -Flow('depthVelModel',['kirchhoffSection','vrms'], +Flow('depthVelModel','vrms', + ''' + time2depth velocity=$SOURCES nz=%i dz=%g intime=y twoway=y | + put label1=Depth unit1=Km + ''' % (nz,dz)) + +# Depth section +Flow('depthSection',['kirchhoffSection','vrms'], ''' time2depth velocity=${SOURCES[1]} nz=%i dz=%g intime=y | put label1=Depth unit1=Km ''' % (nz,dz)) +# Iterative time to depth conversion +Flow('t2d tt2d xt2d ft2d gt2d ct2d','depthVelModel vrms', + ''' + tdconvert niter=3 verb=y cgiter=2000 eps=5 shape=y rect1=4 rect2=15 dix=${SOURCES[1]} + t0=${TARGETS[1]} x0=${TARGETS[2]} f0=${TARGETS[3]} grad=${TARGETS[4]} cost=${TARGETS[5]} + ''') + # Use aliases to split building Alias('model','multiLayerDataCube.rsf') # Show error message if fail - message = ''' SCosncript Failed build ''' From 0f3d075cc565d916a37af879cc853085aebfb655 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sun, 7 Jun 2020 20:34:39 -0300 Subject: [PATCH 07/30] Study about model update #5 --- experiments/modelUpdateStudy/README.md | 1 + experiments/modelUpdateStudy/SConstruct | 83 +++++++++++++++++++++ experiments/modelUpdateStudy/errorReport.py | 54 ++++++++++++++ 3 files changed, 138 insertions(+) create mode 100644 experiments/modelUpdateStudy/README.md create mode 100644 experiments/modelUpdateStudy/SConstruct create mode 100644 experiments/modelUpdateStudy/errorReport.py diff --git a/experiments/modelUpdateStudy/README.md b/experiments/modelUpdateStudy/README.md new file mode 100644 index 0000000..ce30460 --- /dev/null +++ b/experiments/modelUpdateStudy/README.md @@ -0,0 +1 @@ +#### Study about model update and interface interpolation with cubic spline diff --git a/experiments/modelUpdateStudy/SConstruct b/experiments/modelUpdateStudy/SConstruct new file mode 100644 index 0000000..b050921 --- /dev/null +++ b/experiments/modelUpdateStudy/SConstruct @@ -0,0 +1,83 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# SConstruct (Madagascar Script) +# +# Purpose: Study of interface interpolation with cubic spline +# and velocity model update. +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 07/06/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Madagascar package +from rsf.proj import * + +# Error report function +import atexit +from errorReport import print_build_failures + +from random import uniform, seed + +seed() +xmax = 6.0 +zmax = 2.0 + +def arr2str(array,sep=' '): + ''' + Convert a tuple into a comma separeted string + ''' + return string.join(map(str,array),sep) + +# Velocity Model update +for i in range(2): + + interfaces = 'iterfaces-%i' %i + dipsfile = 'dipsfile-%i' %i + modelfile = 'model-%i' %i + ascInterfaces = 'layers-%i.asc' %i + layersFile = 'layers-%i' %i + + layers = [1.65,1.85,1.55,1.65] + layers = map(lambda x: x+uniform(-0.1,0.1),layers) + + velocities = [1.508,2.0] + + vstr=arr2str(velocities,',') + + n1 = 4 #len(layers[0]) + n2 = 1 #len(layers) + + Flow(ascInterfaces,None, + ''' + echo %s + n1=%d n2=%d o1=0 d1=%g + data_format=ascii_float in=$TARGET + ''' % (arr2str(layers),n1,n2,xmax/(n1-1))) + + Flow(layersFile,ascInterfaces,'dd form=native') + + d = 0.0101 # non-round for reproducibility + + Flow(interfaces,layersFile, + 'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d))) + Flow(dipsfile,interfaces,'deriv scale=y') + + Flow(modelfile,interfaces, + ''' + unif2 d1=%g n1=%d v00=%s + ''' % (d,int(1.5+zmax/d),vstr)) + +# Show error message if fail +message = ''' +SCosncript Failed build +''' + +atexit.register(print_build_failures,message) +End() diff --git a/experiments/modelUpdateStudy/errorReport.py b/experiments/modelUpdateStudy/errorReport.py new file mode 100644 index 0000000..8509563 --- /dev/null +++ b/experiments/modelUpdateStudy/errorReport.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# errorReport.py (Madagascar Recipe) +# +# Purpose: SCons error report function. +# +# Important!: It should be called from a SConstruct +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 01/06/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Selfdoc string +''' +Python function to report errors in SCons + +If build fail report specific message. Import this module and add the following lineis to your SConstruct: + +import atexit +from madagascarRecipes.errorReport import print_build_failures + +Call error message by the following: + +atexit.register(print_build_failures,message) + +It will call function print_build_failures to report message variable that should be a string with the error message to be reported +''' + +if __name__ == "__main__": + print(__doc__) + exit() + +__author__="Rodolfo Dirack " +__version__="1.0" + +def print_build_failures(message): + ''' + Print error message when build fail + :param message: string, error message + ''' + from SCons.Script import GetBuildFailures + for bf in GetBuildFailures(): + print("BUILD ERROR REPORT".center(40,'*')) + print("FAILED BUILD: %s" % (bf.node)) + print("MESSAGE REPORT".center(40,'*')) + print("%s" % (message)) + From b7154f9ac67ee0111ae7549ac5e7789fbce7c2b3 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sat, 8 Aug 2020 18:30:08 -0300 Subject: [PATCH 08/30] Ray tracer study #5 Use sfrays2 to ray tracing on a smooth velocity model. --- experiments/rayTracer/SConstruct | 94 ++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 experiments/rayTracer/SConstruct diff --git a/experiments/rayTracer/SConstruct b/experiments/rayTracer/SConstruct new file mode 100644 index 0000000..bca7594 --- /dev/null +++ b/experiments/rayTracer/SConstruct @@ -0,0 +1,94 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# SConstruct (Madagascar Script) +# +# Purpose: Study of ray tracing. +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 08/08/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Madagascar package +from rsf.proj import * + +raysfile='rays' +modelfile='model' +overlay='ovelay' +escfile='esc' + +xmax = 6.0 +zmax = 2.0 + +def arr2str(array,sep=' '): + ''' + Convert a tuple into a comma separeted string + ''' + return string.join(map(str,array),sep) + +# Velocity Model +interfaces = 'iterfaces' +dipsfile = 'dipsfile' +ascInterfaces = 'layers.asc' +layersFile = 'layers' + +layers = [0.65,0.85,0.55,0.65] +#layers = map(lambda x: x+uniform(-0.1,0.1),layers) +velocities = [1.508,2.0] + +vstr=arr2str(velocities,',') + +n1 = 4 #len(layers[0]) +n2 = 1 #len(layers) + +Flow(ascInterfaces,None, + ''' + echo %s + n1=%d n2=%d o1=0 d1=%g + data_format=ascii_float in=$TARGET + ''' % (arr2str(layers),n1,n2,xmax/(n1-1))) + +Flow(layersFile,ascInterfaces,'dd form=native') + +d = 0.00101 # non-round for reproducibility + +Flow(interfaces,layersFile,'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d))) + +Flow(dipsfile,interfaces,'deriv scale=y') + +Flow(modelfile,interfaces, + ''' + unif2 d1=%g n1=%d v00=%s | + smooth rect1=5 rect2=5 repeat=5 + ''' % (d,int(1.5+zmax/d),vstr)) + + + +# Shoot ray in the model +Flow(raysfile,modelfile, + ''' + rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=n + ''') + +# Keep escape values of rays (z,x,t,angle) +Flow(escfile,modelfile, + ''' + rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=y + ''') + +# plot the model +Plot(modelfile,'grey color=j scalebar=y label1=Depth unit1=km label2=Position unit2=km barlabel=Velocity barunit=km/s barreverse=y title=Model allpos=y') + +# plot the ray +Plot(raysfile,'graph transp=y yreverse=y min1=0 max1=2 min2=0 max2=6 wantaxis=n wanttitle=n scalebar=y plotcol=7 plotfat=3') + +# overlay model and ray +Result(overlay,[modelfile, raysfile],'Overlay') + +End() From 138946cd018608a1a4ce9a70d6e9bc97bdfa3aba Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sat, 8 Aug 2020 18:31:57 -0300 Subject: [PATCH 09/30] Add error report library to model update SConstruct #5 --- experiments/modelUpdateStudy/SConstruct | 35 ++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/experiments/modelUpdateStudy/SConstruct b/experiments/modelUpdateStudy/SConstruct index b050921..4aaace6 100644 --- a/experiments/modelUpdateStudy/SConstruct +++ b/experiments/modelUpdateStudy/SConstruct @@ -36,15 +36,18 @@ def arr2str(array,sep=' '): return string.join(map(str,array),sep) # Velocity Model update -for i in range(2): +for i in range(1): interfaces = 'iterfaces-%i' %i dipsfile = 'dipsfile-%i' %i modelfile = 'model-%i' %i ascInterfaces = 'layers-%i.asc' %i layersFile = 'layers-%i' %i + raysfile = 'rays-%i' %i + huygens = 'huygens-%i' %i + overlay = 'overlay-%i' %i - layers = [1.65,1.85,1.55,1.65] + layers = [0.65,0.85,0.55,0.65] layers = map(lambda x: x+uniform(-0.1,0.1),layers) velocities = [1.508,2.0] @@ -63,7 +66,7 @@ for i in range(2): Flow(layersFile,ascInterfaces,'dd form=native') - d = 0.0101 # non-round for reproducibility + d = 0.00101 # non-round for reproducibility Flow(interfaces,layersFile, 'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d))) @@ -71,9 +74,33 @@ for i in range(2): Flow(modelfile,interfaces, ''' - unif2 d1=%g n1=%d v00=%s + unif2 d1=%g n1=%d v00=%s ''' % (d,int(1.5+zmax/d),vstr)) + Flow(raysfile,modelfile, + ''' + rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=n + ''') + +# Flow(huygens,modelfile, +# ''' +# smooth repeat=3 rect1=5 rect2=3 | +# hwt2d xsou=0 zsou=3 nt=1000 d1=0.01 ng=1801 og=90 dg=1 rays=y +# ''') + + # plot the model + Plot(modelfile,'grey color=j scalebar=y label1=Depth unit1=km label2=Position unit2=km barlabel=Velocity barunit=km/s barreverse=y title=Model allpos=y') + + # plot the ray + Plot(raysfile,'graph transp=y yreverse=y min1=0 max1=2 min2=0 max2=6 wantaxis=n wanttitle=n scalebar=y plotcol=7 plotfat=3') + + # plot the ray +# Plot(huygens,'graph transp=n yreverse=y min1=0 max1=2 min2=0 max2=6 wantaxis=n wanttitle=n scalebar=y plotcol=7 plotfat=3') + + + # overlay model and ray + Result(overlay,[modelfile, raysfile],'Overlay') + # Show error message if fail message = ''' SCosncript Failed build From e1db12230da0defc76f13ffda1834c51fa5caae8 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sat, 8 Aug 2020 20:45:31 -0300 Subject: [PATCH 10/30] Simulate a point source #5 Someone can simulated a point source and store rays escape values in a 'txt' file with sfdisfil: ~$ sfdisfil col=4 < esc.rsf > out.txt The collumns of escape values file are z,x,t,angle coordinates for each ray. --- experiments/rayTracer/SConstruct | 13 ++++++++----- experiments/rayTracer/out.txt | 10 ++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) create mode 100644 experiments/rayTracer/out.txt diff --git a/experiments/rayTracer/SConstruct b/experiments/rayTracer/SConstruct index bca7594..a4bb3d3 100644 --- a/experiments/rayTracer/SConstruct +++ b/experiments/rayTracer/SConstruct @@ -22,6 +22,9 @@ raysfile='rays' modelfile='model' overlay='ovelay' escfile='esc' +a0=-45 +nr=10 +amax=45 xmax = 6.0 zmax = 2.0 @@ -38,7 +41,7 @@ dipsfile = 'dipsfile' ascInterfaces = 'layers.asc' layersFile = 'layers' -layers = [0.65,0.85,0.55,0.65] +layers = [0.65,0.65,0.65,0.65] #layers = map(lambda x: x+uniform(-0.1,0.1),layers) velocities = [1.508,2.0] @@ -73,14 +76,14 @@ Flow(modelfile,interfaces, # Shoot ray in the model Flow(raysfile,modelfile, ''' - rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=n - ''') + rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=%g amax=%g nr=%i escvar=n + ''' % (a0,amax,nr)) # Keep escape values of rays (z,x,t,angle) Flow(escfile,modelfile, ''' - rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=y - ''') + rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=%g amax=%g nr=%i escvar=y + ''' % (a0,amax,nr)) # plot the model Plot(modelfile,'grey color=j scalebar=y label1=Depth unit1=km label2=Position unit2=km barlabel=Velocity barunit=km/s barreverse=y title=Model allpos=y') diff --git a/experiments/rayTracer/out.txt b/experiments/rayTracer/out.txt new file mode 100644 index 0000000..e972a1f --- /dev/null +++ b/experiments/rayTracer/out.txt @@ -0,0 +1,10 @@ + 0: -0.005877 2.234 0.76 -32.42 + 4: -0.00739 2.449 0.69 -24.96 + 8: -0.01286 2.624 0.65 -17.79 + 12: -0.00511 2.778 0.62 -11.11 + 16: -0.008527 2.925 0.61 -3.878 + 20: -0.008527 3.075 0.61 3.878 + 24: -0.00511 3.222 0.62 11.11 + 28: -0.01286 3.376 0.65 17.79 + 32: -0.00739 3.551 0.69 24.96 + 36: -0.005877 3.766 0.76 32.42 From 8706183b758a48c5e724252947badcc38c52c32e Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 12 Aug 2020 11:45:11 -0300 Subject: [PATCH 11/30] Original velocity continuation recipes #5 --- .../velocityInversion/velconRecipes/uncert.py | 163 ++++++++++++++++ .../velocityInversion/velconRecipes/velcon.py | 176 ++++++++++++++++++ 2 files changed, 339 insertions(+) create mode 100644 experiments/velocityInversion/velconRecipes/uncert.py create mode 100644 experiments/velocityInversion/velconRecipes/velcon.py diff --git a/experiments/velocityInversion/velconRecipes/uncert.py b/experiments/velocityInversion/velconRecipes/uncert.py new file mode 100644 index 0000000..0dddbe5 --- /dev/null +++ b/experiments/velocityInversion/velconRecipes/uncert.py @@ -0,0 +1,163 @@ +from rsf.proj import * +import velcon + +def uncert(data, # data name + nv, # continuation steps + v0, # initial velocity + dv, # velocity step + nx, # lateral dimension + nh, # number of offsets + padt, # time padding + padt2, # extra time padding + padx=None, # lateral padding + v1=None, # other velocity + n1=None, # time extent + dt=0.004, # time sampling + dx=0.01, # midpoint sampling + units='km', # lateral units + vslope=None, # semblance muting + vx0=0, # semblance muting + x0=0, # lateral origin + rect1=10, # vertical smoothing + rect2=10): # lateral smoothing + + velcon.velcon(data,nv,v0,dv,nx,nh,padt,padt2,padx,v1,n1,dt,dx,units,vslope,vx0,x0,rect1,rect2) + + vlf=data+'-vlf' + vlf2=data+'-vlf1' + npk = data+'-npk' + + # To estimate uncertainty: measure dt/dv, measure dv, multiply + ref = data+'-ref' + Flow(ref,[vlf,npk], + ''' + transp | refer ref=${SOURCES[1]} | window n1=9 f1=%d | transp + ''' % (nv-5)) + + if vslope: + refer = ''' + mutter x0=%g v0=%g half=n | + transp | refer ref=${SOURCES[1]} | transp + ''' % (vx0,vslope) + else: + refer = ''' + transp | refer ref=${SOURCES[1]} | transp + ''' + + Flow(ref+'2',[vlf2,npk],refer) + + dtdv = data+'-dtdv' + Flow(dtdv,ref, + ''' + dip n4=0 rect1=%d rect2=5 rect3=%d | + window n2=1 min2=0 | + scale dscale=%g + ''' % (rect1,rect2,dt/dv)) + Flow(dtdv+'0',[dtdv,npk], + 'math vel=${SOURCES[1]} output="vel*x1*input*input" ') + Result(dtdv, + ''' + grey title="Structural Sensitivity in T" + label1=Time unit1=s label2="Lateral Position" unit2=%s + scalebar=y color=j allpos=y barlabel="dt/dv" barunit="s\^\s75 2\s100 \_/%s" + ''' % (units,units)) + + dxdv = data+'-dxdv' + Flow(dxdv,ref, + ''' + transp plane=13 | + dip n4=0 rect1=%d rect2=5 rect3=%d | + window n2=1 min2=0 | + scale dscale=%g | transp + ''' % (rect1,rect2,dx/dv)) + Result(dxdv, + ''' + grey title="Structural Sensitivity in X" + label1=Time unit1=s label2="Lateral Position" unit2=%s + scalebar=y color=j barlabel="dx/dv" barunit=s + ''' % units) + + ddv = data+'-ddv' + Flow(ref+'3',ref+'2','stack norm=n') + Flow(ddv,[ref+'2',ref+'3'], + ''' + math output="x2*x2*input" | + stack norm=n | + add mode=d ${SOURCES[1]} | + math output="sqrt(input)" + ''') + Result(ddv, + ''' + grey title="Velocity Uncertainty" allpos=y + label1=Time unit1=s label2="Lateral Position" unit2=%s + scalebar=y color=j barlabel=Velocity barunit="%s/s" + ''' % (units,units)) + + unc = data+'-unc' + Flow(unc,[dtdv,ddv],'math dv=${SOURCES[1]} output="abs(0.5*input*dv)" ') + Result(unc, + ''' + grey title="Structural Uncertainty" + scalebar=y color=j allpos=y + label1=Time unit1=s label2="Lateral Position" unit2=%s + barlabel="Vertical Uncertainty (s)" + ''' % units) + Flow(unc+'2',[dxdv,ddv], + 'math dv=${SOURCES[1]} output="abs(0.5*input*dv)" ') + Result(unc+'2', + ''' + grey title="Structural Uncertainty" + scalebar=y color=j allpos=y + label1=Time unit1=s label2="Lateral Position" unit2=%s + barlabel="Lateral Uncertainty (%s)" + ''' % (units,units)) + + ddip=data+'-ddip' + Flow(ddip,data, + ''' + window | + dip rect1=%d rect2=%d rect3=1 + ''' % (rect1,rect2)) + dpwd=data+'-dpwd' + Flow(dpwd,[data,ddip],'window | pwd dip=${SOURCES[1]}') + + data2 = data+'2' + Flow(data2,dpwd,'window n4=1 | spray axis=2 n=1') + + mig2=data2+'-mig' + Flow(mig2,data2,'preconstkirch vel=%g' % v0) + + cip2=data2+'-cip' + Flow(cip2,mig2,velcon.mig2cip) + + if padx: + pad2=data2+'-pad' + Flow(pad2,cip2,'pad n3=%d' % padx) + else: + pad2=cip2 + + ckx2=data2+'-ckx' + vlf2=data2+'-vlf' + Flow(ckx2,pad2,'cosft sign3=1') + Flow(vlf2,ckx2, + ''' + fourvc nv=%d dv=%g v0=%g pad=%d pad2=%d | + cosft sign3=-1 | window n3=%d + ''' % (nv,dv,v0,padt,padt2,nx)) + + + vlf22=data2+'-vlf2' + Flow(vlf22,pad2, + ''' + transp plane=23 memsize=500 | + fourvc2 nv=%d dv=%g v0=%g pad=%d pad2=%d | + window n2=%d | transp plane=23 memsize=500 + ''' % (nv,dv,v0,padt,padt2,nx)) + + + foc=data2+'-foc' + Flow(foc,[vlf2,vlf22], + ''' + focus rect1=%d rect3=%d | + math vlf=${SOURCES[1]} output=vlf/input + ''' % (2*rect1,2*rect2)) diff --git a/experiments/velocityInversion/velconRecipes/velcon.py b/experiments/velocityInversion/velconRecipes/velcon.py new file mode 100644 index 0000000..2296e80 --- /dev/null +++ b/experiments/velocityInversion/velconRecipes/velcon.py @@ -0,0 +1,176 @@ +from rsf.proj import * +import dix + +mig2cip = None + +def velcon(data, # data name + nv, # continuation steps + v0, # initial velocity + dv, # velocity step + nx, # lateral dimension + nh, # number of offsets + padt, # time padding + padt2, # extra time padding + padx=None, # lateral padding + v1=None, # other velocity + n1=None, # time extent + dt=0.004, # time sampling + dx=0.01, # midpoint sampling + units='km', # lateral units + vslope=None, # semblance muting + vx0=0, # semblance muting + x0=0, # lateral origin + srect1=3, # semblance vertical smoothing + srect2=1, # semblance lateral smoothing + rect1=10, # vertical smoothing + rect2=10): # lateral smoothing + '''Velocity continuation''' + global mig2cip + + vm = v0+0.5*nv*dv + + mig=data+'-mig' + Flow(mig,data, + ''' + halfint inv=y adj=y | + preconstkirch vel=%g + ''' % v0,split=[4,nh]) + + if n1: + mig2cip = ''' + transp plane=34 memsize=500 | + transp plane=23 | window n1=%d + ''' % n1 + else: + mig2cip = ''' + transp plane=34 memsize=500 | + transp plane=23 + ''' + n1=100 + + cip=data+'-cip' + Flow(cip,mig,mig2cip) + + if padx: + pad=data+'-pad' + Flow(pad,cip,'pad n3=%d' % padx) + else: + pad=cip + padx=nx + + ckx=data+'-ckx' + vlf=data+'-vlf' + + ckx2=data+'-ckx2' + vlf2=data+'-vlf2' + + Flow(ckx,pad,'cosft sign3=1 | put o4=0') + Flow(ckx+'v',ckx, + ''' + fourvc nv=%d dv=%g v0=%g pad=%d pad2=%d verb=y + ''' % (nv,dv,v0,padt,padt2), + split=[3,padx]) + Flow(vlf,ckx+'v', + ''' + cosft sign3=-1 | window n3=%d + ''' % nx) + + Flow(ckx2,pad, + ''' + halfint inv=y adj=y | + math output="input*input" | + halfint adj=y | + cosft sign3=1 | put o4=0 + ''') + Flow(ckx2+'v',ckx2, + ''' + fourvc nv=%d dv=%g v0=%g pad=%d pad2=%d verb=y + ''' % (nv,dv,v0,padt,padt2), + split=[3,padx]) + Flow(vlf2,ckx2+'v', + ''' + cosft sign3=-1 | window n3=%d | clip2 lower=0 + ''' % nx) + + sem = data+'-sem' + Flow(sem,[vlf,vlf2], + ''' + mul $SOURCE | + divn den=${SOURCES[1]} rect1=%d rect3=%d + ''' % (srect1,srect2)) + + vlf1 = data+'-vlf1' + + Flow(vlf1,pad, + ''' + transp plane=23 memsize=1000 | + fourvc2 nv=%d dv=%g v0=%g pad=%d pad2=%d | + window n2=%d | transp plane=23 memsize=1000 + ''' % (nv,dv,v0,padt,padt2,nx)) + + if v1: + Flow(mig+'1',data,'preconstkirch vel=%g' % v1,split=[4,nh]) + Flow(cip+'1',mig+'1',mig2cip) + + migr = data+'-migr' + + Flow(migr,cip,'stack norm=y') + Plot(migr,'grey title=Migration0') + + Flow(migr+'1',cip+'1','stack norm=y') + Plot(migr+'1','grey title=Migration1') + + vlfr = data+'-vlfr' + + Flow(vlfr,vlf,'window n2=1 min2=%g' % v1) + Plot(vlfr,'grey title="Velocity Continuation 0 -> 1" ') + + Result(migr,[migr,migr+'1',vlfr],'SideBySideAniso') + + + if vslope: + pick = ''' + mutter x0=%g v0=%g half=n | + scale axis=2 | pick rect1=%d rect2=%d + ''' % (vx0,vslope,rect1,rect2) + else: + pick = ''' + scale axis=2 | pick rect1=%d rect2=%d + ''' % (rect1,rect2) + + npk = data+'-npk' + Flow(npk,sem,pick) + Plot(npk, + ''' + grey pclip=100 color=j bias=%g + scalebar=y title="Picked Migration Velocity" + label1=Time unit1=s label2="Lateral Position" unit2=%s + barlabel=Velocity barunit="%s/s" barreverse=y + ''' % (vm,units,units)) + + fmg = data+'-fmg' + Flow(fmg,[vlf,npk],'slice pick=${SOURCES[1]}') + + Result(fmg, + ''' + grey title=Slice label1=Time unit1=s + label2="Lateral Position" unit2=%s + ''' % units) + + agc = data+'-agc' + Flow(agc,fmg,'agc rect1=200') + Plot(fmg,agc, + ''' + grey title="Time-Migrated Image" label1="Time" + unit1=s label2="Lateral Position" unit2=%s + ''' % units) + Result(agc, + ''' + grey title=Picked pclip=98 label1="Time" + unit1=s label2="Lateral Position" unit2=%s + ''' % units) + + Result(fmg+'2',[fmg,npk],'SideBySideAniso',vppen='txscale=1.2') + + Flow(agc+'2',[sem,npk],'slice pick=${SOURCES[1]} | window') + From 7cf63ec902c80202c6d5b3675176690b9fa66502 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 12 Aug 2020 11:45:58 -0300 Subject: [PATCH 12/30] Pre-stack velocity analisys with velocity continuation #5 Use program sffourvc to do velocity analisys by velocity continuation. --- experiments/velocityInversion/SConscript2 | 65 +++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 experiments/velocityInversion/SConscript2 diff --git a/experiments/velocityInversion/SConscript2 b/experiments/velocityInversion/SConscript2 new file mode 100644 index 0000000..6460519 --- /dev/null +++ b/experiments/velocityInversion/SConscript2 @@ -0,0 +1,65 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# SConscript2 (Madagascar Script) +# +# Purpose: Study about velocity continuation. +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 11/08/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Madagascar package +from rsf.proj import * + +# Error report function +import atexit +from madagascarRecipes.errorReport import print_build_failures + +# CRE recipe +from madagascarRecipes.pefInterpolation import pefInterpolation as pefin +from madagascarRecipes.kimodel import multiLayerModel as mlmod +from madagascarRecipes.kimodel import kirchoffNewtonModeling as kinewmod +from madagascarRecipes.velocityAnalisys import velocityAnalisys as nmoStack + +# Velocity continuation recipe +from rsf.recipes.velcon import velcon + +Flow('dataTransposed','multiLayerDataCube','transp plane=23 | transp plane=34') + +velcon(data='dataTransposed', # data name + nv=100, # continuation steps + v0=1.5, # initial velocity + dv=0.01, # velocity step + nx=201, # lateral dimension + nh=161, # number of offsets + padt=1024, # time padding + padt2=2048, # extra time padding + padx=None, # lateral padding + v1=None, # other velocity + n1=None, # time extent + dt=0.004, # time sampling + dx=0.025, # midpoint sampling + units='km', # lateral units + vslope=None, # semblance muting + vx0=0, # semblance muting + x0=0, # lateral origin + srect1=3, # semblance vertical smoothing + srect2=1, # semblance lateral smoothing + rect1=10, # vertical smoothing + rect2=10) # lateral smoothing + + +# Show error message if fail +message = ''' +SCosncript Failed build +''' + +atexit.register(print_build_failures,message) +End() From 99d933a36eeb162f0123cf9d14e6adeff04ae760 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 19 Aug 2020 14:47:25 -0300 Subject: [PATCH 13/30] Study about velocity continuation Resolve #5 Usage example of post stack velocity continuation with sfvczo --- experiments/velocityInversion/SConscript2 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/experiments/velocityInversion/SConscript2 b/experiments/velocityInversion/SConscript2 index 6460519..03aefd1 100644 --- a/experiments/velocityInversion/SConscript2 +++ b/experiments/velocityInversion/SConscript2 @@ -33,6 +33,7 @@ from rsf.recipes.velcon import velcon Flow('dataTransposed','multiLayerDataCube','transp plane=23 | transp plane=34') +# Pre-stack velocity continuation velcon(data='dataTransposed', # data name nv=100, # continuation steps v0=1.5, # initial velocity @@ -55,6 +56,15 @@ velcon(data='dataTransposed', # data name rect1=10, # vertical smoothing rect2=10) # lateral smoothing +# Post-stack velocity continuation +Flow('velcon','stackedSection', + ''' + pad beg2=200 end2=200 | cosft sign2=1 | + stolt vel=1.5 | + vczo v0=1.5 dv=0.01 nv=100 verb=y | + transp plane=23 | cosft sign2=-1 | + window min2=0 max2=5 + ''') # Show error message if fail message = ''' From 69725157900bdbc6fdb48fd6b4bab45d9511cc04 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 19 Aug 2020 14:49:37 -0300 Subject: [PATCH 14/30] Original files of velocity continuation programs --- experiments/velocityContinuation/Mfourvc.c | 246 ++++++++++++++++++++ experiments/velocityContinuation/Mfourvc2.c | 245 +++++++++++++++++++ experiments/velocityContinuation/Mvczo.c | 172 ++++++++++++++ experiments/velocityContinuation/README.md | 1 + 4 files changed, 664 insertions(+) create mode 100644 experiments/velocityContinuation/Mfourvc.c create mode 100644 experiments/velocityContinuation/Mfourvc2.c create mode 100644 experiments/velocityContinuation/Mvczo.c create mode 100644 experiments/velocityContinuation/README.md diff --git a/experiments/velocityContinuation/Mfourvc.c b/experiments/velocityContinuation/Mfourvc.c new file mode 100644 index 0000000..b5b3f23 --- /dev/null +++ b/experiments/velocityContinuation/Mfourvc.c @@ -0,0 +1,246 @@ +/* Prestack velocity continuation. */ +/* + Copyright (C) 2004 University of Texas at Austin + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include "fint1.h" +#include +#include + +#ifdef SF_HAS_FFTW +#include +#endif + +int main(int argc, char* argv[]) +{ + fint1 str, istr; + bool verb; + int i1,i2, n1,n2,n3, nw, nx,ny,nv,nh, ix,iy,iv,ih, next; + float d1,o1,d2,o2, eps, w,x,y,k, v0,v2,v,v1,dv, dx,dy, h0,dh,h, t, x0, y0, dw; + float *trace, *strace; + sf_complex *ctrace, *ctrace2, **cstack, shift; + char *time, *space, *unit; + size_t len; + sf_file in, out; + +#ifdef SF_HAS_FFTW + fftwf_plan forw, invs; +#else + kiss_fftr_cfg forw, invs; +#endif + + sf_init (argc,argv); + in = sf_input("in"); + out = sf_output("out"); + + if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input"); + if (!sf_histint(in,"n2",&nh)) sf_error("No n2= in input"); + if (!sf_histint(in,"n3",&nx)) sf_error("No n3= in input"); + if (!sf_histint(in,"n4",&ny)) ny=1; + + if (!sf_getfloat("eps",&eps)) eps=0.01; /* regularization */ + if (!sf_getint("pad",&n2)) n2=n1; /* padding for stretch */ + if (!sf_getint("pad2",&n3)) n3=2*kiss_fft_next_fast_size((n2+1)/2); + /* padding for FFT */ + + if (!sf_getbool("verb",&verb)) verb=false; + /* verbosity flag */ + + nw = n3/2+1; + + if (!sf_histfloat(in,"o1",&o1)) o1=0.; + o2 = o1*o1; + + if(!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input"); + d2 = o1+(n1-1)*d1; + d2 = (d2*d2 - o2)/(n2-1); + + if (!sf_getint("nv",&nv)) sf_error("Need nv="); + /* velocity steps */ + if (!sf_getfloat("dv",&dv)) sf_error("Need dv="); + /* velocity step size */ + if (!sf_getfloat("v0",&v0) && + !sf_histfloat(in,"v0",&v0)) sf_error("Need v0="); + /*( v0 starting velocity )*/ + + if(!sf_histfloat(in,"o2",&h0)) sf_error("No o2= in input"); + if(!sf_histfloat(in,"d2",&dh)) sf_error("No d2= in input"); + + if(!sf_histfloat(in,"d3",&dx)) sf_error("No d3= in input"); + if(!sf_histfloat(in,"d4",&dy)) dy=dx; + if(!sf_histfloat(in,"o3",&x0)) x0=0.; + if(!sf_histfloat(in,"o4",&y0)) y0=0.; + + sf_putfloat(out,"o2",v0+dv); + sf_putfloat(out,"d2",dv); + sf_putint(out,"n2",nv); + + sf_putstring(out,"label2","Velocity"); + + if (NULL != (time = sf_histstring(in,"unit1")) && + NULL != (space = sf_histstring(in,"unit2"))) { + len = strlen(time)+strlen(space)+2; + unit = sf_charalloc(len); + snprintf(unit,len,"%s/%s",space,time); + sf_putstring(out,"unit2",unit); + free(time); + free(space); + } + + dx *= 2.*SF_PI; + dy *= 2.*SF_PI; + x0 *= 2.*SF_PI; + y0 *= 2.*SF_PI; + + dw = 16*SF_PI/(d2*n3); /* 2pi * 8 */ + + trace = sf_floatalloc(n1); + strace = sf_floatalloc(n3); + ctrace = sf_complexalloc(nw); + ctrace2 = sf_complexalloc(nw); + cstack = sf_complexalloc2(nw,nv); + +#ifdef SF_HAS_FFTW + forw = fftwf_plan_dft_r2c_1d(n3, strace, (fftwf_complex *) ctrace, + FFTW_ESTIMATE); + invs = fftwf_plan_dft_c2r_1d(n3, (fftwf_complex *) ctrace2, strace, + FFTW_ESTIMATE); +#else + forw = kiss_fftr_alloc(n3,0,NULL,NULL); + invs = kiss_fftr_alloc(n3,1,NULL,NULL); +#endif + if (NULL == forw || NULL == invs) sf_error("FFT allocation error"); + + + if (!sf_getint("extend",&next)) next=4; + /* trace extension */ + str = fint1_init(next,n1,0); + istr = fint1_init(next,n2,0); + + for (iy=0; iy < ny; iy++) { + y = y0+iy*dy; + y *= y; + for (ix=0; ix < nx; ix++) { + x = x0+ix*dx; + x *= x; + + if (verb) sf_warning("wavenumber %d of %d and %d of %d;", + ix+1,nx, iy+1,ny); + k = (x+y) * 0.5; + + for (iv=0; iv < nv; iv++) { + for (i1=0; i1 < nw; i1++) { + cstack[iv][i1] = sf_cmplx(0.,0.); + } + } + + for (ih=0; ih < nh; ih++) { + h = h0 + ih*dh; + h *= h * 0.5; + + sf_floatread(trace,n1,in); + for (i1=0; i1 < n1; i1++) { + trace[i1] /= (n1*nh); + } + fint1_set(str,trace); + + for (i2=0; i2 < n2; i2++) { + t = o2+i2*d2; + t = sqrtf(t); + t = (t-o1)/d1; + i1 = t; + if (i1 >= 0 && i1 < n1) { + strace[i2] = fint1_apply(str,i1,t-i1,false); + } else { + strace[i2] = 0.; + } + } + + for (i2=n2; i2 < n3; i2++) { + strace[i2] = 0.; + } + +#ifdef SF_HAS_FFTW + fftwf_execute(forw); +#else + kiss_fftr(forw,strace, (kiss_fft_cpx *) ctrace); +#endif + + for (iv=0; iv < nv; iv++) { + v = v0 + (iv+1)*dv; + v1 = h * (1./(v*v) - 1./(v0*v0)); + v2 = k * ((v0*v0) - (v*v)); + + for (i2=1; i2 < nw; i2++) { + w = i2*dw; + w = v2/w+(v1-0.125*o2)*w; + shift = sf_cmplx(cosf(w),sinf(w)); + +#ifdef SF_HAS_COMPLEX_H + cstack[iv][i2] += ctrace[i2] * shift; +#else + cstack[iv][i2] = sf_cadd(cstack[iv][i2], + sf_cmul(ctrace[i2],shift)); +#endif + } /* w */ + } /* v */ + } /* h */ + + for (iv=0; iv < nv; iv++) { + + ctrace2[0] = sf_cmplx(0.0f,0.0f); /* dc */ + + for (i2=1; i2 < nw; i2++) { + w = i2*dw; + w *= 0.125*o2; + shift = sf_cmplx(cosf(w),sinf(w)); + +#ifdef SF_HAS_COMPLEX_H + ctrace2[i2] = cstack[iv][i2] * shift; +#else + ctrace2[i2] = sf_cmul(cstack[iv][i2],shift); +#endif + } + +#ifdef SF_HAS_FFTW + fftwf_execute(invs); +#else + kiss_fftri(invs,(const kiss_fft_cpx *) ctrace2, strace); +#endif + + fint1_set(istr,strace); + + for (i1=0; i1 < n1; i1++) { + t = o1+i1*d1; + t = t*t; + t = (t-o2)/d2; + i2 = t; + if (i2 >= 0 && i2 < n2) { + trace[i1] = fint1_apply(istr,i2,t-i2,false); + } else { + trace[i1] = 0.; + } + } + + sf_floatwrite (trace,n1,out); + } /* v 2 */ + } /* x */ + } /* y */ + if (verb) sf_warning("."); + + exit (0); +} diff --git a/experiments/velocityContinuation/Mfourvc2.c b/experiments/velocityContinuation/Mfourvc2.c new file mode 100644 index 0000000..d1db8e9 --- /dev/null +++ b/experiments/velocityContinuation/Mfourvc2.c @@ -0,0 +1,245 @@ +/* Velocity continuation with semblance computation. */ +/* + Copyright (C) 2004 University of Texas at Austin + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include +#include +#include "fint1.h" + +int main(int argc, char* argv[]) +{ + fint1 str, istr; + int i1,i2, n1,n2,n3, ix,iv,ih, ib, ie, nb, nx,nv,nh, nw, next; + float d1,o1,d2,o2, eps, w,x,k, v0,v2,v,v1,dv, dx, h0,dh,h, num, den, t, dw; + float *trace=NULL, *strace=NULL, ***stack=NULL, ***stack2=NULL, ***cont=NULL, **image=NULL; + sf_complex *ctrace=NULL, *ctrace0=NULL, shift; + char *time=NULL, *space=NULL, *unit=NULL; + size_t len; + static kiss_fftr_cfg forw, invs; + sf_file in=NULL, out=NULL; + bool sembl; + + sf_init (argc,argv); + in = sf_input("in"); + out = sf_output("out"); + + if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input"); + if (!sf_histint(in,"n2",&nx)) sf_error("No n2= in input"); + if (!sf_histint(in,"n3",&nh)) sf_error("No n3= in input"); + + if (!sf_getint("nb",&nb)) nb=2; + if (!sf_getfloat("eps",&eps)) eps=0.01; + if (!sf_getint("pad",&n2)) n2=n1; + if (!sf_getint("pad2",&n3)) n3=2*kiss_fft_next_fast_size((n2+1)/2); + + nw = n3/2+1; + forw = kiss_fftr_alloc(n3,0,NULL,NULL); + invs = kiss_fftr_alloc(n3,1,NULL,NULL); + if (NULL == forw || NULL == invs) + sf_error("KISS FFT allocation error"); + + if (!sf_histfloat(in,"o1",&o1)) o1=0.; + o2 = o1*o1; + + if(!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input"); + d2 = o1+(n1-1)*d1; + d2 = (d2*d2 - o2)/(n2-1); + + if (!sf_getint("nv",&nv)) sf_error("Need nv="); + if (!sf_getfloat("dv",&dv)) sf_error("Need dv="); + if (!sf_getfloat("v0",&v0) && + !sf_histfloat(in,"v0",&v0)) sf_error("Need v0="); + + if (!sf_getbool("semblance",&sembl)) sembl=true; + /* if y, compute semblance; if n, stack */ + + if(!sf_histfloat(in,"o3",&h0)) sf_error("No o2= in input"); + if(!sf_histfloat(in,"d3",&dh)) sf_error("No d2= in input"); + if(!sf_histfloat(in,"d2",&dx)) sf_error("No d3= in input"); + + sf_putfloat(out,"o3",v0+dv); + sf_putfloat(out,"d3",dv); + sf_putint(out,"n3",nv); + + sf_putstring(out,"label3","Velocity"); + + if (NULL != (time = sf_histstring(in,"label1")) && + NULL != (space = sf_histstring(in,"label2"))) { + len = strlen(time)+strlen(space)+2; + unit = sf_charalloc(len); + snprintf(unit,len,"%s/%s",space,time); + sf_putstring(out,"unit3",unit); + free(time); + free(space); + } + + dx = 2*SF_PI/(2*kiss_fft_next_fast_size(nx-1)*dx); + dw = 16*SF_PI/(d2*n3); /* 2pi * 8 */ + + stack = sf_floatalloc3(n1,nx,nv); + stack2 = sf_floatalloc3(n1,nx,nv); + cont = sf_floatalloc3(n1,nx,nv); + image = sf_floatalloc2(n1,nx); + trace = sf_floatalloc(n1); + strace = sf_floatalloc(n3); + ctrace = sf_complexalloc(nw); + ctrace0 = sf_complexalloc(nw); + + if (!sf_getint("extend",&next)) next=4; + /* trace extension */ + str = fint1_init(next,n1,0); + istr = fint1_init(next,n2,0); + + for (i1=0; i1 < n1*nx*nv; i1++) { + stack[0][0][i1] = 0.; + stack2[0][0][i1] = 0.; + } + + sf_cosft_init(nx); + + for (ih=0; ih < nh; ih++) { + sf_warning("offset %d of %d;",ih+1,nh); + + h = h0 + ih*dh; + h *= h * 0.5; + + sf_floatread(image[0],n1*nx,in); + + for (i1=0; i1 < n1; i1++) { + sf_cosft_frw(image[0],i1,n1); + } + + for (ix=0; ix < nx; ix++) { + x = ix*dx; + x *= x; + + k = x * 0.5; + + fint1_set(str,image[ix]); + + for (i2=0; i2 < n2; i2++) { + t = o2+i2*d2; + t = sqrtf(t); + t = (t-o1)/d1; + i1 = t; + if (i1 >= 0 && i1 < n1) { + strace[i2] = fint1_apply(str,i1,t-i1,false); + } else { + strace[i2] = 0.; + } + } + + for (i2=n2; i2 < n3; i2++) { + strace[i2] = 0.; + } + + kiss_fftr(forw,strace, (kiss_fft_cpx *) ctrace0); + + for (iv=0; iv < nv; iv++) { + v = v0 + (iv+1)* dv; + + v1 = h * (1./(v*v) - 1./(v0*v0)); + v2 = k * ((v0*v0) - (v*v)); + + ctrace[0]=sf_cmplx(0.,0.); /* dc */ + + for (i2=1; i2 < nw; i2++) { + w = i2*dw; + w = v2/w+(v1-0.125*o2)*w; + shift = sf_cmplx(cosf(w),sinf(w)); + +#ifdef SF_HAS_COMPLEX_H + ctrace[i2] = ctrace0[i2] * shift; +#else + ctrace[i2] = sf_cmul(ctrace0[i2],shift); +#endif + } /* w */ + + kiss_fftri(invs,(const kiss_fft_cpx *) ctrace, strace); + + fint1_set(istr,strace); + + for (i1=0; i1 < n1; i1++) { + t = o1+i1*d1; + t = t*t; + t = (t-o2)/d2; + i2 = t; + if (i2 >= 0 && i2 < n2) { + cont[iv][ix][i1] = fint1_apply(istr,i2,t-i2,false); + } else { + cont[iv][ix][i1] = 0.; + } + } + } /* v */ + } /* x */ + + for (iv=0; iv < nv; iv++) { + for (i1=0; i1 < n1; i1++) { + sf_cosft_inv(cont[0][0],i1+iv*nx*n1,n1); + } + } + + for (iv=0; iv < nv; iv++) { + for (ix=0; ix < nx; ix++) { + for (i1=0; i1 < n1; i1++) { + t = cont[iv][ix][i1]; + stack[iv][ix][i1] += t; + stack2[iv][ix][i1] += t*t; + } /* i1 */ + } /* x */ + } /* v */ + } /* h */ + sf_warning("."); + + for (iv=0; iv < nv; iv++) { + for (ix=0; ix < nx; ix++) { + for (i1=0; i1 < n1; i1++) { + ib = i1-nb > 0? i1-nb: 0; + ie = i1+nb+1 < n1? i1+nb+1: n1; + + num = 0.; + den = 0.; + + if (sembl) { + + for (i2=ib; i2 < ie; i2++) { + t = stack[iv][ix][i2]; + num += t*t; + den += stack2[iv][ix][i2]; + } + + den *= nh; + + trace[i1] = den > 0.? num/den: 0.; + } else { + + for (i2=ib; i2 < ie; i2++) { + t = stack[iv][ix][i2]; + num += t; + } + + den = nh; + trace[i1] = num/(den+ FLT_EPSILON); + } + } + sf_floatwrite (trace,n1,out); + } + } + + exit(0); +} diff --git a/experiments/velocityContinuation/Mvczo.c b/experiments/velocityContinuation/Mvczo.c new file mode 100644 index 0000000..b336172 --- /dev/null +++ b/experiments/velocityContinuation/Mvczo.c @@ -0,0 +1,172 @@ +/* Post-stack 2-D velocity continuation. */ +/* + Copyright (C) 2004 University of Texas at Austin + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ +#include + +#ifdef SF_HAS_FFTW +#include +#endif + +int main(int argc, char* argv[]) +{ + sf_map4 str; + bool verb; + int i1,i2, n1,n2,n3, nw, nx,nv, ix,iv, iy,ny; + float d1,o1,d2,o2, eps, w,x, v0,v2,v,dv, dx, t, x0, dw; + float *trace, *strace, *t2; + sf_complex *ctrace, *ctrace2, shift; + sf_file in, out; + +#ifdef SF_HAS_FFTW + fftwf_plan forw, invs; +#else + kiss_fftr_cfg forw, invs; +#endif + + sf_init (argc,argv); + in = sf_input("in"); + out = sf_output("out"); + + if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input"); + if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input"); + if (!sf_histint(in,"n2",&nx)) sf_error("No n2= in input"); + ny = sf_leftsize(in,2); + + if (!sf_getfloat("eps",&eps)) eps=0.01; /* regularization */ + if (!sf_getint("pad",&n2)) n2=n1; /* padding for stretch */ + if (!sf_getint("pad2",&n3)) n3=2*kiss_fft_next_fast_size((n2+1)/2); + /* padding for FFT */ + + if (!sf_getbool("verb",&verb)) verb=true; + /* verbosity flag */ + + nw = n3/2+1; + + strace = sf_floatalloc(n3); + ctrace = sf_complexalloc(nw); + ctrace2 = sf_complexalloc(nw); + +#ifdef SF_HAS_FFTW + forw = fftwf_plan_dft_r2c_1d(n3, strace, (fftwf_complex *) ctrace, + FFTW_ESTIMATE); + invs = fftwf_plan_dft_c2r_1d(n3, (fftwf_complex *) ctrace2, strace, + FFTW_ESTIMATE); +#else + forw = kiss_fftr_alloc(n3,0,NULL,NULL); + invs = kiss_fftr_alloc(n3,1,NULL,NULL); +#endif + if (NULL == forw || NULL == invs) sf_error("FFT allocation error"); + + if (!sf_histfloat(in,"o1",&o1)) o1=0.; + o2 = o1*o1; + + if(!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input"); + d2 = o1+(n1-1)*d1; + d2 = (d2*d2 - o2)/(n2-1); + dw = 16*SF_PI/(d2*n3); /* 2pi * 8 */ + + if (!sf_getint("nv",&nv)) sf_error("Need nv="); + /* velocity steps */ + if (!sf_getfloat("dv",&dv)) sf_error("Need dv="); + /* velocity step size */ + if (!sf_getfloat("v0",&v0) && + !sf_histfloat(in,"v0",&v0)) sf_error("Need v0="); + /*( v0 starting velocity )*/ + + if(!sf_histfloat(in,"d2",&dx)) sf_error("No d2= in input"); + if(!sf_histfloat(in,"o2",&x0)) x0=0.; + + sf_putfloat(out,"o2",v0+dv); + sf_putfloat(out,"d2",dv); + sf_putint(out,"n2",nv); + + sf_putstring(out,"label2","Velocity"); + + sf_shiftdim(in, out, 2); + + dx *= 2.*SF_PI; + x0 *= 2.*SF_PI; + + trace = sf_floatalloc(n1); + t2 = sf_floatalloc(n2); + + str = sf_stretch4_init (n1, o1, d1, n2, eps); + + for (i2=0; i2 < n2; i2++) { + t = o2+i2*d2; + t2[i2] = sqrtf(t); + } + + sf_stretch4_define (str,t2); + + for (iy=0; iy < ny; iy++) { + for (ix=0; ix < nx; ix++) { + if (verb) sf_warning("wavenumber %d of %d (%d);", ix+1,nx,iy); + + x = x0+ix*dx; + x *= x; + x *= 0.5; + + sf_floatread(trace,n1,in); + for (i1=0; i1 < n1; i1++) { + trace[i1] /= n1; + } + sf_stretch4_invert (false,str,strace,trace); + for (i2=n2; i2 < n3; i2++) { + strace[i2] = 0.; + } + +#ifdef SF_HAS_FFTW + fftwf_execute(forw); +#else + kiss_fftr(forw,strace, (kiss_fft_cpx *) ctrace); +#endif + + for (iv=0; iv < nv; iv++) { + v = v0 + (iv+1)*dv; + v2 = x * ((v0*v0) - (v*v)); + + ctrace2[0] = sf_cmplx(0.0f,0.0f); /* dc */ + + for (i2=1; i2 < nw; i2++) { + w = i2*dw; + w = v2/w; + + shift = sf_cmplx(cosf(w),sinf(w)); + +#ifdef SF_HAS_COMPLEX_H + ctrace2[i2] = ctrace[i2] * shift; +#else + ctrace2[i2] = sf_cmul(ctrace[i2],shift); +#endif + } /* w */ + +#ifdef SF_HAS_FFTW + fftwf_execute(invs); +#else + kiss_fftri(invs,(const kiss_fft_cpx *) ctrace2, strace); +#endif + sf_stretch4_apply(false,str,strace,trace); + sf_floatwrite (trace,n1,out); + } /* v */ + } /* x */ + if (verb) sf_warning("."); + } /* y */ + + exit (0); +} diff --git a/experiments/velocityContinuation/README.md b/experiments/velocityContinuation/README.md new file mode 100644 index 0000000..4de558a --- /dev/null +++ b/experiments/velocityContinuation/README.md @@ -0,0 +1 @@ +#### Original source codes of velocity continuation programs From b7eccc34de5071049d3410685b830c0ca37ba328 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 19 Aug 2020 15:08:23 -0300 Subject: [PATCH 15/30] Scratch of velocity inversion algorithm #2 In this directory scratch a velocity inversion algorithm using velocity continuation and diffraction simulation --- experiments/algorithmVelocityInversion/README.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 experiments/algorithmVelocityInversion/README.md diff --git a/experiments/algorithmVelocityInversion/README.md b/experiments/algorithmVelocityInversion/README.md new file mode 100644 index 0000000..bec8360 --- /dev/null +++ b/experiments/algorithmVelocityInversion/README.md @@ -0,0 +1 @@ +#### Scratch of velocity inversion algorithm From 111b06cd094af540d2d5a0e8a3365915f1ad4d85 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 19 Aug 2020 17:11:20 -0300 Subject: [PATCH 16/30] SConstruct of velocity inversion algorithm #2 It uses MadagascarRecipes installed in local Madagascar package to be called from SCons. Install recipes source files in $RSFSRC/book/Recipes directory and install it in your local Madagascar version with the command 'sudo scons install' called from $RSFSRC. --- .../algorithmVelocityInversion/SConstruct | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 experiments/algorithmVelocityInversion/SConstruct diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct new file mode 100644 index 0000000..882d1e9 --- /dev/null +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -0,0 +1,91 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# SConstruct (Madagascar Script) +# +# Purpose: Scratch of velocity inversion algorithm. +# +# Site: http://www.dirackslounge.online +# +# Version 1.0 +# +# Programer: Rodolfo A. C. Neves (Dirack) 19/08/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +# Madagascar package +from rsf.proj import * + +# Personal Madagascar recipes +from rsf.recipes.kimodel import multiLayerModel as mlmod +from rsf.recipes.kimodel import kirchoffNewtonModeling as kinewmod +from rsf.recipes.velocityAnalisys import velocityAnalisys as nmoStack + +# Velocity continuation recipe +from rsf.recipes.velcon import velcon + +xmax = 6.0 +zmax = 2.0 + +layers = ((0.30,0.50,0.20,0.30), + (1.65,1.85,1.55,1.65)) + +velocities = (1.508, + 1.690, + 2.0) + +# Generate multi layer model and data cube +mlmod(interfaces='interfaces', + dipsfile='interfacesDip', + modelfile='mod1', + xmax=xmax, + zmax=zmax, + layers=layers, + velocities=velocities) + +kinewmod(reflectors='interfaces', + reflectorsDip='interfacesDip', + filename='multiLayerDataCube', + velocities=velocities) + +Flow('dataTransposed','multiLayerDataCube','transp plane=23 | transp plane=34') + +# Pre-stack velocity continuation +velcon(data='dataTransposed', # data name + nv=100, # continuation steps + v0=1.5, # initial velocity + dv=0.01, # velocity step + nx=201, # lateral dimension + nh=161, # number of offsets + padt=1024, # time padding + padt2=2048, # extra time padding + padx=None, # lateral padding + v1=None, # other velocity + n1=None, # time extent + dt=0.004, # time sampling + dx=0.025, # midpoint sampling + units='km', # lateral units + vslope=None, # semblance muting + vx0=0, # semblance muting + x0=0, # lateral origin + srect1=3, # semblance vertical smoothing + srect2=1, # semblance lateral smoothing + rect1=10, # vertical smoothing + rect2=10) # lateral smoothing + +# Post-stack velocity continuation +Flow('velocityCube','stackedSection', + ''' + pad beg2=200 end2=200 | cosft sign2=1 | + stolt vel=1.5 | + vczo v0=1.5 dv=0.01 nv=100 verb=y | + transp plane=23 | cosft sign2=-1 | + window min2=0 max2=5 + ''') + +# Use aliases to split building +Alias('model','multiLayerDataCube.rsf') + +End() From 07711139af862cc2045a1a9266c997671db4dfc2 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Wed, 19 Aug 2020 18:31:20 -0300 Subject: [PATCH 17/30] Correct parameters of modeling recipe #2 Some parameters and function names have changed in MadagascarRecipe package. --- .../algorithmVelocityInversion/SConstruct | 31 ++++++++++++++++--- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 882d1e9..2a9d9d7 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -20,8 +20,8 @@ from rsf.proj import * # Personal Madagascar recipes from rsf.recipes.kimodel import multiLayerModel as mlmod -from rsf.recipes.kimodel import kirchoffNewtonModeling as kinewmod -from rsf.recipes.velocityAnalisys import velocityAnalisys as nmoStack +from rsf.recipes.kimodel import kirchhoffNewtonModeling as kinewmod +from rsf.recipes.velocityAnalysis import velocityAnalysis as nmoStack # Velocity continuation recipe from rsf.recipes.velcon import velcon @@ -48,7 +48,13 @@ mlmod(interfaces='interfaces', kinewmod(reflectors='interfaces', reflectorsDip='interfacesDip', filename='multiLayerDataCube', - velocities=velocities) + velocities=velocities, + nt=1001, + dt=0.004, + ns=201, + ds=0.025, + nh=161, + dh=0.025) Flow('dataTransposed','multiLayerDataCube','transp plane=23 | transp plane=34') @@ -74,7 +80,21 @@ velcon(data='dataTransposed', # data name srect2=1, # semblance lateral smoothing rect1=10, # vertical smoothing rect2=10) # lateral smoothing - + +# velocity analysis, NMO correction and stack +nmoStack(dataCube='multiLayerDataCube', + pick='vnmo', + stack='stackedSection', + vrms='vrms', + v0=1.5, + dv=0.01, + nv=100, + vel0=1.5, + rect1=15, + rect2=40, + rect3=3, + dt=0.004) + # Post-stack velocity continuation Flow('velocityCube','stackedSection', ''' @@ -86,6 +106,9 @@ Flow('velocityCube','stackedSection', ''') # Use aliases to split building +# Call it 'scons alias' to build +# target identified by alias Alias('model','multiLayerDataCube.rsf') +Alias('stack','velocityCube.rsf') End() From e44af5f58a45876a5d85aba85874261f969c01c6 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 14:28:28 -0300 Subject: [PATCH 18/30] Iterative picking #2 Iterative picking of the stacked section with sfipick --- experiments/algorithmVelocityInversion/SConstruct | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 2a9d9d7..5c5da66 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -105,6 +105,12 @@ Flow('velocityCube','stackedSection', window min2=0 max2=5 ''') +# Reflector iterative Picking +Flow('reflectorPickedPoints.txt','stackedSection', + ''' + ipick | tr '\t' ' ' + ''') + # Use aliases to split building # Call it 'scons alias' to build # target identified by alias From f5bfe9f8440336475e2c540bb537e3f99da4badd Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 14:49:57 -0300 Subject: [PATCH 19/30] Scratch of t0 and m0 files building #2 T0 and m0 files store (t0,m0) picked pairs from stacked section. Those points will be used to simulate diffraction hyperbolas in the stacked section with centre in (t0,m0) coordinates. --- .../algorithmVelocityInversion/SConstruct | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 5c5da66..a7bcd0e 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -108,7 +108,21 @@ Flow('velocityCube','stackedSection', # Reflector iterative Picking Flow('reflectorPickedPoints.txt','stackedSection', ''' - ipick | tr '\t' ' ' + ipick + ''') + +# Build t0 coordinates file +Flow('t0s.txt','reflectorPickedPoints.txt', + ''' + tr '\\t' ' ' | + /usr/bin/cut -d' ' -f1 + ''') + +# Build m0 coordinates file +Flow('m0s.txt','reflectorPickedPoints.txt', + ''' + tr '\\t' ' ' | + /usr/bin/cut -d' ' -f2 ''') # Use aliases to split building From 8e3d72bff8ca575b21351d5fca5697a3574908a5 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 15:03:30 -0300 Subject: [PATCH 20/30] Velocities file #2 Build the velocities file to store layers velocity that will be used to build the simulated diffraction hyperbolas --- experiments/algorithmVelocityInversion/SConstruct | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index a7bcd0e..63de4bd 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -125,6 +125,16 @@ Flow('m0s.txt','reflectorPickedPoints.txt', /usr/bin/cut -d' ' -f2 ''') +# Build velocity files +Flow('v.asc',None, + ''' + echo %s %s + n1=2 o1=0 d1=1 + data_format=ascii_float in=$TARGET + ''' % (velocities[0],velocities[1])) + +Flow('v','v.asc','dd form=native') + # Use aliases to split building # Call it 'scons alias' to build # target identified by alias From 22b4edab4fed0143fb7788f716009d704721fad0 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 21:08:16 -0300 Subject: [PATCH 21/30] Shell Script to format picked points to ascii #2 The Shell Script 'ascFormat.sh' formats sfipick picked points to ascii format readable by sfdd that converts it to RSF format. It reads points from a sfipick generated file and converts it to RSF. If user pass 1 to the script, it cuts the first collumn of the sfipick output and builds a RSF file with t0's coordinates. If user pass 2 to the script, it cuts the second collumn of the sfipick output and builds a RSF file with m0's coordinates. --- .../algorithmVelocityInversion/ascFormat.sh | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100755 experiments/algorithmVelocityInversion/ascFormat.sh diff --git a/experiments/algorithmVelocityInversion/ascFormat.sh b/experiments/algorithmVelocityInversion/ascFormat.sh new file mode 100755 index 0000000..3cee240 --- /dev/null +++ b/experiments/algorithmVelocityInversion/ascFormat.sh @@ -0,0 +1,33 @@ +#!/bin/bash +# +# ascFormat.sh (Shell Script) +# +# Purpose: Convert sfipick picked points +# to Madagascar readable ascii format. +# Important! $1=1 cut and format the first collumn of the file (t0 time) +# If $1=2 cut and format the second collumn of the file (m0 CMP) +# +# Site: https://dirack.github.io +# +# Version 1.0 +# +# Programer: Rodolfo A C Neves (Dirack) 20/08/2020 +# +# Email: rodolfo_profissional@hotmail.com +# +# License: GPL-3.0 . + +INPUT=$(cat "-" | tr '\t' ' ' | cut -d' ' -f"$1") +NUMBER_OF_POINTS=$(echo "$INPUT" | wc -l) +FORMATED_INPUT=$(echo "$INPUT" | tr '\n' ' ') + +FILE="t0s" +if [ "$1" == "2" ] +then + FILE="m0s" +else + echo "Error: $(basename $0): \$1 should be equal to 1 or 2!" + exit 1 +fi + +echo "${FORMATED_INPUT}n1=$NUMBER_OF_POINTS d1=1 o1=0 data_format=ascii_float in=${FILE}.asc" From c397c3ecc929e690e88971a21e50fb00b4a9c067 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 21:13:07 -0300 Subject: [PATCH 22/30] Modify SConstruct to use ascFormat.sh script #2 Use the ascFormat.sh Shell Script to convert txt file generated by sfipick to t0s and m0s ascii files readable by sfdd program that converts ascii to RSF. --- .../algorithmVelocityInversion/SConstruct | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 63de4bd..77862d3 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -111,18 +111,18 @@ Flow('reflectorPickedPoints.txt','stackedSection', ipick ''') -# Build t0 coordinates file -Flow('t0s.txt','reflectorPickedPoints.txt', +# Next step is done with 'ascFormat.sh' Shell Script +# please check if the script have permissions to execute +# Build t0 coordinates file (pass 1 to generate t0s file) +Flow('t0s.asc','reflectorPickedPoints.txt', ''' - tr '\\t' ' ' | - /usr/bin/cut -d' ' -f1 + ./ascFormat.sh 1 ''') -# Build m0 coordinates file -Flow('m0s.txt','reflectorPickedPoints.txt', +# Build m0 coordinates file (pass 2 to generate m0s file) +Flow('m0s.asc','reflectorPickedPoints.txt', ''' - tr '\\t' ' ' | - /usr/bin/cut -d' ' -f2 + ./ascFormat.sh 2 ''') # Build velocity files From bd35f714f64430fcc358d396dbed23fa731b9f72 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Thu, 20 Aug 2020 21:18:44 -0300 Subject: [PATCH 23/30] Convert t0 and m0 ascii files to RSF #2 Convert files with t0s and m0s coordinates from ascii to RSF --- experiments/algorithmVelocityInversion/SConstruct | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 77862d3..89cd5f2 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -119,12 +119,16 @@ Flow('t0s.asc','reflectorPickedPoints.txt', ./ascFormat.sh 1 ''') +Flow('t0s','t0s.asc','sfdd form=native') + # Build m0 coordinates file (pass 2 to generate m0s file) Flow('m0s.asc','reflectorPickedPoints.txt', ''' ./ascFormat.sh 2 ''') +Flow('m0s','m0s.asc','sfdd form=native') + # Build velocity files Flow('v.asc',None, ''' From d64da486824e1fa18160acd3ddf2ccf6ba213a6a Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Fri, 21 Aug 2020 16:45:41 -0300 Subject: [PATCH 24/30] Correct bug in parameters input of ascFormat.sh #2 Correct if to check if user passed 1 or 2 from command line in order to decide which collumn to cut from input file to format --- experiments/algorithmVelocityInversion/ascFormat.sh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/experiments/algorithmVelocityInversion/ascFormat.sh b/experiments/algorithmVelocityInversion/ascFormat.sh index 3cee240..607dad6 100755 --- a/experiments/algorithmVelocityInversion/ascFormat.sh +++ b/experiments/algorithmVelocityInversion/ascFormat.sh @@ -21,8 +21,10 @@ INPUT=$(cat "-" | tr '\t' ' ' | cut -d' ' -f"$1") NUMBER_OF_POINTS=$(echo "$INPUT" | wc -l) FORMATED_INPUT=$(echo "$INPUT" | tr '\n' ' ') -FILE="t0s" -if [ "$1" == "2" ] +if [ "$1" == "1" ] +then + FILE="t0s" +elif [ "$1" == "2" ] then FILE="m0s" else From 7df188e6d9c0e4b0e0869816f56c823d08a23f86 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Fri, 21 Aug 2020 16:47:36 -0300 Subject: [PATCH 25/30] Diffraction simulation hyperbolas above one reflector #2 Use program sfdiffsim to simulated hyperbolas over picked pairs in t0s.rsf and m0s.rsf files with a given velocity --- experiments/algorithmVelocityInversion/SConstruct | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 89cd5f2..f61f94f 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -129,15 +129,11 @@ Flow('m0s.asc','reflectorPickedPoints.txt', Flow('m0s','m0s.asc','sfdd form=native') -# Build velocity files -Flow('v.asc',None, +# Diffraction simulation in stacked section +Flow(['returnedSection','diffSection'],['stackedSection','t0s','m0s'], ''' - echo %s %s - n1=2 o1=0 d1=1 - data_format=ascii_float in=$TARGET - ''' % (velocities[0],velocities[1])) - -Flow('v','v.asc','dd form=native') + diffsim diff=${TARGETS[1]} t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y + ''' % (velocities[0])) # Use aliases to split building # Call it 'scons alias' to build From b797477861bea2289529341f63dc04ecbf325fba Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Fri, 21 Aug 2020 17:47:01 -0300 Subject: [PATCH 26/30] Shell Script receives ascii filename #2 It is easy to pass the filename to a Shell Script instead of pass it through stdin --- experiments/algorithmVelocityInversion/ascFormat.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/experiments/algorithmVelocityInversion/ascFormat.sh b/experiments/algorithmVelocityInversion/ascFormat.sh index 607dad6..e7c5172 100755 --- a/experiments/algorithmVelocityInversion/ascFormat.sh +++ b/experiments/algorithmVelocityInversion/ascFormat.sh @@ -32,4 +32,4 @@ else exit 1 fi -echo "${FORMATED_INPUT}n1=$NUMBER_OF_POINTS d1=1 o1=0 data_format=ascii_float in=${FILE}.asc" +echo "${FORMATED_INPUT}n1=$NUMBER_OF_POINTS d1=1 o1=0 data_format=ascii_float in=$2" From 0fca09cbb88ab554fcf60b270d4d344c55af90b5 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Fri, 21 Aug 2020 17:51:46 -0300 Subject: [PATCH 27/30] Generate simulated diffraction hyperbolas in a loop #2 Generate a set of diffraction hyperbolas in a loop. Each iteration is assigned to a reflector defined in layers vector --- .../algorithmVelocityInversion/SConstruct | 74 +++++++++++-------- 1 file changed, 45 insertions(+), 29 deletions(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index f61f94f..3ee561a 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -105,35 +105,51 @@ Flow('velocityCube','stackedSection', window min2=0 max2=5 ''') -# Reflector iterative Picking -Flow('reflectorPickedPoints.txt','stackedSection', - ''' - ipick - ''') - -# Next step is done with 'ascFormat.sh' Shell Script -# please check if the script have permissions to execute -# Build t0 coordinates file (pass 1 to generate t0s file) -Flow('t0s.asc','reflectorPickedPoints.txt', - ''' - ./ascFormat.sh 1 - ''') - -Flow('t0s','t0s.asc','sfdd form=native') - -# Build m0 coordinates file (pass 2 to generate m0s file) -Flow('m0s.asc','reflectorPickedPoints.txt', - ''' - ./ascFormat.sh 2 - ''') - -Flow('m0s','m0s.asc','sfdd form=native') - -# Diffraction simulation in stacked section -Flow(['returnedSection','diffSection'],['stackedSection','t0s','m0s'], - ''' - diffsim diff=${TARGETS[1]} t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y - ''' % (velocities[0])) +# Loop over reflectors +numberOfReflectors = len(layers) +section = 'stackedSection' + +for i in range(numberOfReflectors): + + reflectorPickedPoints = 'reflectorPickedPoints-%i.txt' % i + t0sFile = 't0s-%i' % i + t0sAscii = 't0s-%i.asc' %i + m0sFile = 'm0s-%i' % i + m0sAscii = 'm0s-%i.asc' % i + returnedSection = 'returnedSection-%i' % i + diffSection = 'diffSection-%i' % i + + # Reflector iterative Picking + Flow(reflectorPickedPoints,section, + ''' + ipick + ''') + + # Next step is done with 'ascFormat.sh' Shell Script + # please check if the script have permissions to execute + # Build t0 coordinates file (pass 1 to generate t0s file) + Flow(t0sAscii,reflectorPickedPoints, + ''' + ./ascFormat.sh 1 %s + ''' % (t0sAscii)) + + Flow(t0sFile,t0sAscii,'sfdd form=native') + + # Build m0 coordinates file (pass 2 to generate m0s file) + Flow(m0sAscii,reflectorPickedPoints, + ''' + ./ascFormat.sh 2 %s + ''' % (m0sAscii)) + + Flow(m0sFile,m0sAscii,'sfdd form=native') + + # Diffraction simulation in stacked section + Flow([returnedSection,diffSection],[section,t0sFile,m0sFile], + ''' + diffsim diff=${TARGETS[1]} t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y + ''' % (velocities[i])) + + section = returnedSection # Use aliases to split building # Call it 'scons alias' to build From 122c5cdb8f583f10bbcc6bd3c6ccc1f20edb22f9 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sat, 22 Aug 2020 12:39:28 -0300 Subject: [PATCH 28/30] Add simulated sections #2 Add diffraction simulated sections and migrate it --- .../algorithmVelocityInversion/SConstruct | 33 ++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 3ee561a..c7c1ed2 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -22,6 +22,7 @@ from rsf.proj import * from rsf.recipes.kimodel import multiLayerModel as mlmod from rsf.recipes.kimodel import kirchhoffNewtonModeling as kinewmod from rsf.recipes.velocityAnalysis import velocityAnalysis as nmoStack +from rsf.recipes.diff import diffmig as diff # Velocity continuation recipe from rsf.recipes.velcon import velcon @@ -146,11 +147,41 @@ for i in range(numberOfReflectors): # Diffraction simulation in stacked section Flow([returnedSection,diffSection],[section,t0sFile,m0sFile], ''' - diffsim diff=${TARGETS[1]} t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y + diffsim diff=${TARGETS[1]} aperture=0.5 + t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y ''' % (velocities[i])) section = returnedSection +Flow('diffAdd',['diffSection-0','diffSection-1'],'sfadd ${SOURCES[1]} scale=1,1') + +diff("diff", + 'diffAdd', + v0=1.5, + nv=20, + dv=0.1, + nx=201, + padx=1000, + nt=1000, + tmin=0, + tmax=4, + rect1=10, + rect2=10, + srect1=1, + srect2=3, + vslope=None, + units='Km', + f1=1, + j3=1, + dx=0.025, + x0=0, + beg1=0, + frect1=0, + frect2=0, + an=1, + nout=2048, + vx0=None) + # Use aliases to split building # Call it 'scons alias' to build # target identified by alias From cf2ceecffaffaad0b3196550a11082ca8b51aa14 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sun, 23 Aug 2020 16:18:01 -0300 Subject: [PATCH 29/30] Increase hyperbolas aperture to 1Km #2 --- experiments/algorithmVelocityInversion/SConstruct | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index c7c1ed2..00d42a5 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -147,7 +147,7 @@ for i in range(numberOfReflectors): # Diffraction simulation in stacked section Flow([returnedSection,diffSection],[section,t0sFile,m0sFile], ''' - diffsim diff=${TARGETS[1]} aperture=0.5 + diffsim diff=${TARGETS[1]} aperture=1 t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y ''' % (velocities[i])) From 9e71670c30c905b8524bdbfda403717d99cc3fb2 Mon Sep 17 00:00:00 2001 From: Rodolfo Dirack Date: Sun, 23 Aug 2020 16:20:04 -0300 Subject: [PATCH 30/30] Increase velocity window of diffraction migration #2 Change velocity sampling and increase number of steps to improve velocity searching for a wider window and better resolution. --- experiments/algorithmVelocityInversion/SConstruct | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/experiments/algorithmVelocityInversion/SConstruct b/experiments/algorithmVelocityInversion/SConstruct index 00d42a5..f9a4b4f 100644 --- a/experiments/algorithmVelocityInversion/SConstruct +++ b/experiments/algorithmVelocityInversion/SConstruct @@ -157,9 +157,9 @@ Flow('diffAdd',['diffSection-0','diffSection-1'],'sfadd ${SOURCES[1]} scale=1,1' diff("diff", 'diffAdd', - v0=1.5, - nv=20, - dv=0.1, + v0=1.4, + nv=100, + dv=0.01, nx=201, padx=1000, nt=1000,