{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# notebook for create init and true test model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import math\n", "\n", "# grid\n", "#R_earth = 6378.1370\n", "R_earth = 6371.0\n", "\n", "rr1=6361 \n", "rr2=6381\n", "tt1=(38.0-0.3)/180*math.pi\n", "tt2=(42.0+0.3)/180*math.pi\n", "pp1=(23.0-0.3)/180*math.pi\n", "pp2=(27.0+0.3)/180*math.pi\n", "\n", "n_rtp = [10,50,50]\n", "n_rtp.reverse()\n", "dr = (rr2-rr1)/(n_rtp[2]-1)\n", "dt = (tt2-tt1)/(n_rtp[1]-1)\n", "dp = (pp2-pp1)/(n_rtp[0]-1)\n", "rr = np.array([rr1 + x*dr for x in range(n_rtp[2])])\n", "tt = np.array([tt1 + x*dt for x in range(n_rtp[1])])\n", "pp = np.array([pp1 + x*dp for x in range(n_rtp[0])])\n", "\n", "# initial model\n", "gamma = 0.0\n", "s0 = 1.0/6.0\n", "slow_p=0.06\n", "ani_p=0.04\n", "\n", "eta_init = np.zeros(n_rtp)\n", "xi_init = np.zeros(n_rtp)\n", "zeta_init = np.zeros(n_rtp)\n", "fun_init = np.zeros(n_rtp)\n", "vel_init = np.zeros(n_rtp)\n", "a_init = np.zeros(n_rtp)\n", "b_init = np.zeros(n_rtp)\n", "c_init = np.zeros(n_rtp)\n", "f_init = np.zeros(n_rtp)\n", "\n", "# true model\n", "eta_true = np.zeros(n_rtp)\n", "xi_true = np.zeros(n_rtp)\n", "zeta_true = np.zeros(n_rtp)\n", "fun_true = np.zeros(n_rtp)\n", "vel_true = np.zeros(n_rtp)\n", "a_true = np.zeros(n_rtp)\n", "b_true = np.zeros(n_rtp)\n", "c_true = np.zeros(n_rtp)\n", "f_true = np.zeros(n_rtp)\n", "\n", "c=0\n", "for ir in range(n_rtp[2]):\n", " for it in range(n_rtp[1]):\n", " for ip in range(n_rtp[0]):\n", " #eta_init[ip,it,ir] = 0.0\n", " #xi_init[ip,it,ir] = 0.0\n", " zeta_init[ip,it,ir] = gamma*math.sqrt(eta_init[ip,it,ir]**2 + xi_init[ip,it,ir]**2)\n", " fun_init[ip,it,ir] = s0\n", " vel_init[ip,it,ir] = 1.0/s0\n", " a_init[ip,it,ir] = 1.0 + 2.0*zeta_init[ip,it,ir]\n", " b_init[ip,it,ir] = 1.0 - 2.0*xi_init[ip,it,ir]\n", " c_init[ip,it,ir] = 1.0 + 2.0*xi_init[ip,it,ir]\n", " f_init[ip,it,ir] = -2.0 * eta_init[ip,it,ir]\n", "\n", " # true model\n", " if (tt[it] >= 38.0/180.0*math.pi and tt[it] <= 42.0/180.0*math.pi \\\n", " and pp[ip] >= 23.0/180.0*math.pi and pp[ip] <= 27.0/180.0*math.pi):\n", " c+=1\n", " sigma = math.sin(2.0*math.pi*(tt[it]-38.0/180.0*math.pi)/(4.0/180.0*math.pi))*math.sin(2.0*math.pi*(pp[ip]-23.0/180.0*math.pi)/(4.0/180.0*math.pi))\n", " else:\n", " sigma = 0.0\n", "\n", " if sigma < 0:\n", " psi = 60.0/180.0*math.pi\n", " elif sigma > 0:\n", " psi = 120.0/180.0*math.pi\n", " else:\n", " psi = 0.0\n", "\n", " eta_true[ip,it,ir] = ani_p*abs(sigma)*math.sin(2.0*psi)\n", " xi_true[ip,it,ir] = ani_p*abs(sigma)*math.cos(2.0*psi)\n", " zeta_true[ip,it,ir] = gamma*math.sqrt(eta_true[ip,it,ir]**2 + xi_true[ip,it,ir]**2)\n", " fun_true[ip,it,ir] = s0/(1.0+sigma*slow_p)\n", " vel_true[ip,it,ir] = 1.0/fun_true[ip,it,ir] \n", " a_true[ip,it,ir] = 1.0 + 2.0*zeta_true[ip,it,ir]\n", " b_true[ip,it,ir] = 1.0 - 2.0*xi_true[ip,it,ir]\n", " c_true[ip,it,ir] = 1.0 + 2.0*xi_true[ip,it,ir]\n", " f_true[ip,it,ir] = -2.0 * eta_true[ip,it,ir]\n", "\n", "\n", "\n", "#r_earth = 6378.1370\n", "print(\"depminmax {} {}\".format(R_earth-rr1,R_earth-rr2))\n", "print(c)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# write out\n", "import h5py\n", "\n", "fout_init = h5py.File('test_model_init.h5', 'w')\n", "fout_true = h5py.File('test_model_true.h5', 'w')\n", "\n", "# write out the arrays eta_init, xi_init, zeta_init, fun_init, a_init, b_init, c_init, f_init\n", "fout_init.create_dataset('eta', data=eta_init.T)\n", "fout_init.create_dataset('xi', data=xi_init.T)\n", "fout_init.create_dataset('zeta', data=zeta_init.T)\n", "fout_init.create_dataset('fun', data=fun_init.T)\n", "fout_init.create_dataset('fac_a', data=a_init.T)\n", "fout_init.create_dataset('fac_b', data=b_init.T)\n", "fout_init.create_dataset('fac_c', data=c_init.T)\n", "fout_init.create_dataset('fac_f', data=f_init.T)\n", "fout_init.create_dataset('vel', data=vel_init.T)\n", "\n", "# writeout the arrays eta_true, xi_true, zeta_true, fun_true, a_true, b_true, c_true, f_true\n", "fout_true.create_dataset('eta', data=eta_true.T)\n", "fout_true.create_dataset('xi', data=xi_true.T)\n", "fout_true.create_dataset('zeta', data=zeta_true.T)\n", "fout_true.create_dataset('fun', data=fun_true.T)\n", "fout_true.create_dataset('fac_a', data=a_true.T)\n", "fout_true.create_dataset('fac_b', data=b_true.T)\n", "fout_true.create_dataset('fac_c', data=c_true.T)\n", "fout_true.create_dataset('fac_f', data=f_true.T)\n", "fout_true.create_dataset('vel', data=vel_true.T)\n", "\n", "fout_init.close()\n", "fout_true.close()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# prepare src station file\n", "\n", "```\n", " 26 1992 1 1 2 43 56.900 1.8000 98.9000 137.00 2.80 8 305644 <- src 　: id_src year month day hour min sec lat lon dep_km mag num_recs id_event\n", " 26 1 PCBI 1.8900 98.9253 1000.0000 P 10.40 18.000 <- arrival : id_src id_rec name_rec lat lon elevation_m phase epicentral_distance_km arrival_time_sec\n", " 26 2 MRPI 1.6125 99.3172 1100.0000 P 50.84 19.400\n", " 26 3 HUTI 2.3153 98.9711 1600.0000 P 57.84 19.200\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import random\n", "random.seed(1145141919810)\n", "\n", "# dummys\n", "year_dummy = 1998\n", "month_dummy = 1\n", "day_dummy = 1\n", "hour_dummy = 0\n", "minute_dummy = 0\n", "second_dummy = 0\n", "mag_dummy = 3.0\n", "id_dummy = 1000\n", "st_name_dummy = 'AAAA'\n", "phase_dummy = 'P'\n", "dist_dummy = 100.0\n", "arriv_t_dummy = 0.0\n", "\n", "tt1deg = tt1 * 180.0/math.pi\n", "tt2deg = tt2 * 180.0/math.pi\n", "pp1deg = pp1 * 180.0/math.pi\n", "pp2deg = pp2 * 180.0/math.pi\n", "\n", "\n", "n_src = 500\n", "n_rec = [30 for x in range(n_src)]\n", "\n", "\n", "lines = []\n", "\n", "nij_src = math.sqrt(n_src)\n", "nij_rec = math.sqrt(n_rec[0])\n", "\n", "pos_src=[]\n", "pos_rec=[]\n", "\n", "\n", "# create receiver coordinates\n", "elev_recs=[]\n", "lon_recs=[]\n", "lat_recs=[]\n", "rec_names=[]\n", "for i in range(n_rec[0]):\n", " #elev_recs.append(random.uniform(-100.0,-100.0)) # elevation in m\n", " #elev_recs.append(0) # elevation in m\n", " #lon_recs .append(random.uniform(pp1deg*1.1,pp2deg*0.9))\n", " #lat_recs .append(random.uniform(tt1deg*1.1,tt2deg*0.9))\n", " rec_names.append(i)\n", " # regularly\n", " elev_recs.append(0.0)\n", " tmp_ilon = i%nij_rec\n", " tmp_ilat = int(i/nij_rec)\n", " lon_recs.append(pp1deg + tmp_ilon*(pp2deg-pp1deg)/nij_rec)\n", " lat_recs.append(tt1deg + tmp_ilat*(tt2deg-tt1deg)/nij_rec)\n", "\n", "\n", "\n", "# create dummy src\n", "for i_src in range(n_src):\n", " # define one point in the domain (rr1 bottom, rr2 top)\n", " # random\n", " #dep = random.uniform((R_earth-rr1)*0.5,(R_earth-rr1)*0.98)\n", " #lon = random.uniform(pp1deg,pp2deg)\n", " #lat = random.uniform(tt1deg,tt2deg)\n", "\n", " # regularl\n", " dep = (R_earth-rr1)*0.9\n", " tmp_ilon = i_src%nij_src\n", " tmp_ilat = int(i_src/nij_src)\n", " lon = pp1deg + tmp_ilon*(pp2deg-pp1deg)/nij_src\n", " lat = tt1deg + tmp_ilat*(tt2deg-tt1deg)/nij_src\n", "\n", " src = [i_src, year_dummy, month_dummy, day_dummy, hour_dummy, minute_dummy, second_dummy, lat, lon, dep, mag_dummy, n_rec[i_src], id_dummy]\n", " lines.append(src)\n", "\n", " pos_src.append([lon,lat,dep])\n", "\n", "\n", " # create dummy station\n", " for i_rec in range(n_rec[i_src]):\n", " #elev_rec = 0.0 #random.uniform(0.0,-10.0) # elevation in m\n", " #lon_rec = random.uniform(pp1deg,pp2deg)\n", " #lat_rec = random.uniform(tt1deg,tt2deg)\n", " # regularly\n", " #elev_rec = -10.0\n", " #tmp_ilon = i_rec%nij_rec\n", " #tmp_ilat = int(i_rec/nij_rec)\n", " #lon_rec = pp1deg + tmp_ilon*(pp2deg-pp1deg)/nij_rec\n", " #lat_rec = tt1deg + tmp_ilat*(tt2deg-tt1deg)/nij_rec\n", "\n", " # \n", " elev_rec = elev_recs[i_rec]\n", " lon_rec = lon_recs[i_rec]\n", " lat_rec = lat_recs[i_rec]\n", " st_name_dummy = rec_names[i_rec]\n", "\n", " rec = [i_src, i_rec, st_name_dummy, lat_rec, lon_rec, elev_rec, phase_dummy, dist_dummy, arriv_t_dummy]\n", " lines.append(rec)\n", "\n", " pos_rec.append([lon_rec,lat_rec,elev_rec])\n", "\n", "\n", "# write out ev_arrivals file\n", "fname = 'src_rec_test.dat'\n", "\n", "with open(fname, 'w') as f:\n", " for line in lines:\n", " for elem in line:\n", " f.write('{} '.format(elem))\n", " f.write('\\n')\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# draw src and rec positions\n", "import matplotlib.pyplot as plt\n", "\n", "for i_src in range(n_src):\n", " plt.scatter(pos_src[i_src][1],pos_src[i_src][0],c='r',marker='o')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot receivers\n", "for i_rec in range(n_rec[0]):\n", " plt.scatter(pos_rec[i_rec][1],pos_rec[i_rec][0],c='b',marker='o')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.1 64-bit ('3.9.1')", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.1" }, "vscode": { "interpreter": { "hash": "02f83e1f4cd9619657a6845405e2dd67c4de23753ba48bca5dce2ebf57b3813a" } } }, "nbformat": 4, "nbformat_minor": 2 }