{
 "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": [
    "# write out in ASCIII\n",
    "\n",
    "#\n",
    "\n",
    "fname_init = 'test_model_init.dat'\n",
    "fname_true = 'test_model_true.dat'\n",
    "\n",
    "\n",
    "# write init model\n",
    "with open(fname_init, 'w') as f:\n",
    "    # write nodes in rtp\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",
    "                # write out eta xi zeta fun fac_a fac_b fac_c fac_f\n",
    "                f.write(\"{}   {}   {}   {}   {}   {}   {}   {}   {}\\n\".format(eta_init[ip,it,ir],xi_init[ip,it,ir],zeta_init[ip,it,ir],fun_init[ip,it,ir],vel_init[ip,it,ir],a_init[ip,it,ir],b_init[ip,it,ir],c_init[ip,it,ir],f_init[ip,it,ir]))\n",
    "\n",
    "\n",
    "# write true model\n",
    "with open(fname_true, 'w') as f:\n",
    "    # write nodes in rtp\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",
    "                # write out eta xi zeta fun fac_a fac_b fac_c fac_f\n",
    "                f.write(\"{}   {}   {}   {}   {}   {}   {}   {}   {}\\n\".format(eta_true[ip,it,ir],xi_true[ip,it,ir],zeta_true[ip,it,ir],fun_true[ip,it,ir],vel_true[ip,it,ir],a_true[ip,it,ir],b_true[ip,it,ir],c_true[ip,it,ir],f_true[ip,it,ir]))\n",
    "\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": "fbd0b2a7df497f398d93ab2f589d8a5daa3108cfb7ff2b90736653cca3aeadc0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}